Statistical analysis of plasma carotenoids, plasma cytokines and immune cells from randomized cross-over USDA inflammation clinical trial. Subjects consumed both low lycopene tomato (yellow) and high lycopene tomato-soy juices (red) for 4 weeks each.
Crossover clinical trial design supplementing individuals with obesity 360 mL of a low carotenoid tomato juice or a high lycopene tomato-soy juice daily. Daily serving of low carotenoid tomato juice consisted of ~1.5mg lycopene/day while high lycopene tomato-soy juice intervention consisted of 54 mg lycopene/day in addition to 210 mg total soy isoflavones/day.
library(tidyverse) # data wrangling
library(readxl) # read in excel files
library(janitor) # clean up names in dataset
library(corrr) # finding correlations
library(rstatix) # stats
library(knitr) # aesthetic table viewing
library(purrr) # create functions
library(kableExtra)
library(ggthemes)
library(ggtext)
library(ggpubr)
library(ggcorrplot) # to make corr plots
library(pheatmap)
# load data
meta_table <- read_excel("CompiledData_Results_Meta.xlsx",
sheet = "metadata_corrected_withsequence")
# clean up variable names
meta_table <- clean_names(meta_table)
str(meta_table)
# convert variables that should be factors to factors
meta_table <- meta_table %>%
mutate(across(.cols = c("patient_id", "period",
"intervention", "intervention_week",
"pre_post", "sex", "sequence"),
.fns = as.factor))
# some stuff came in as characters but should be numeric
meta_table <- meta_table %>%
mutate(across(.cols = c("il_2", "il_10", "il_13", "il_4"),
.fns = as.numeric))
str(meta_table)
# changing factor levels for pre_post
meta_table$pre_post <- factor(meta_table$pre_post,
levels = c("pre", "post"))
levels(meta_table$pre_post)
## [1] "pre" "post"
# Calculate total_cis_lyc, total_lyc, and total_carotenoids
meta_table <- meta_table %>%
rename(n5_cis_lyc = x5_cis_lyc) %>%
mutate(total_cis_lyc = other_cis_lyc + n5_cis_lyc,
total_lyc = all_trans_lyc + total_cis_lyc,
total_carotenoids = lutein + zeaxanthin + b_cryptoxanthin +
a_carotene + b_carotene + total_lyc)
# line plots for each subject at each timepoint
meta_table %>%
ggplot(aes(x = intervention_week, y = all_trans_lyc, color = intervention)) +
geom_point() +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Baseline" = "gray",
"Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_bw() +
labs(x = "Intervention Week",
y = "All-trans-lycopene levels (nmol/L)",
title = "All-trans-lycopene levels in each patient before/after each intervention")
# line plots for each subject at each timepoint
meta_table %>%
ggplot(aes(x = intervention_week, y = total_cis_lyc, color = intervention)) +
geom_point() +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Baseline" = "gray",
"Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id)) +
theme_bw() +
labs(x = "Intervention Week",
y = "Total cis-lycopene levels (nmol/L)",
title = "Total cis-lycopene levels in each patient before/after each intervention")
# line plots for each subject at each timepoint
meta_table %>%
ggplot(aes(x = intervention_week, y = total_lyc, color = intervention)) +
geom_point() +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Baseline" = "gray",
"Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_bw() +
labs(x = "Intervention Week",
y = "Total lycopene levels (nmol/L)",
title = "Total lycopene levels in each patient before/after each intervention")
note Subject 6112 has increased total lycopene levels after yellow intervention. I’ll remove them from statistical analyses.
# create a more specific pre_post_intervention column
meta_table_edited <- meta_table %>%
unite(col = "pre_post_intervention",
c("pre_post","intervention"),
sep = "_",
remove = FALSE)
# make pre_post_intervention column factors
meta_table_edited$pre_post_intervention <- as.factor(meta_table_edited$pre_post_intervention)
# relevel factor columns
meta_table_edited$pre_post_intervention <- factor(meta_table_edited$pre_post_intervention, levels = c("pre_Yellow", "post_Yellow", "pre_Red", "post_Red"))
meta_table_edited$intervention <- factor(meta_table_edited$intervention,
levels = c("Baseline","Yellow", "Red"))
# make legend title
legendtitle_ppintervention <- "Timepoint"
# labels
labs_ppintervention <- c("pre control",
"post control",
"pre Tomato-Soy",
"post Tomato-Soy")
meta_table_edited %>%
filter(intervention != "Baseline") %>%
ggplot(aes(x = intervention, y = total_lyc, fill = pre_post_intervention)) +
geom_boxplot(outlier.shape = NA) +
scale_fill_manual(legendtitle_ppintervention,
values = c("pre_Red" = "#FF9966",
"post_Red" = "#FF3300",
"pre_Yellow" = "#FFFF99",
"post_Yellow" = "yellow"),
labels = labs_ppintervention) +
theme_clean() +
labs(x = "",
y = "Total lycopene levels (nmol/L)",
title = "Total lycopene levels before/after juice interventions")
# ggpubr to make pub ready plot
(total_lyc_levels <- meta_table_edited %>%
filter(intervention != "Baseline") %>%
ggpaired(x = "pre_post", y = "total_lyc", fill = "intervention", line.color = "gray", line.size = 1, facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", "", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", linewidth = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "nmol/L plasma",
title = "Concentration of Lycopene",
subtitle = ""))
# convert meta_table_edited to long format for total lycopene
meta_table_lyc_long <- meta_table_edited %>%
pivot_longer(cols = total_lyc,
names_to = "total_lycopene",
values_to = "nmol_per_L")
lyc_subset <- meta_table_lyc_long %>%
select(patient_id, pre_post_intervention, nmol_per_L) %>%
pivot_wider(names_from = pre_post_intervention,
values_from = nmol_per_L) %>%
mutate(red_FC = post_Red/pre_Red,
yellow_FC = post_Yellow/pre_Yellow)
lyc_subset %>%
summarize(mean_red_FC = mean(red_FC),
mean_yellow_FC = mean(yellow_FC))
## # A tibble: 1 × 2
## mean_red_FC mean_yellow_FC
## <dbl> <dbl>
## 1 2.83 0.966
First, I need to remove subject 6112 from the yellow intervention, since their total lycopene levels is the only one increasing post-Yellow
subj6112_yellow <- meta_table_lyc_long %>%
filter(patient_id == 6112) %>%
filter(intervention == "Yellow")
meta_table_lyc_long_rm_outlier <- anti_join(meta_table_lyc_long,
subj6112_yellow)
Plot histogram
gghistogram(meta_table_lyc_long_rm_outlier$nmol_per_L, bins = 40)
Shapiro’s normality test
# shapiro normality test for total lycopene. for pre yellow, post yellow, pre red and post red
meta_table_lyc_long_rm_outlier %>%
filter(intervention != "Baseline") %>%
group_by(pre_post_intervention) %>%
shapiro_test(vars = "nmol_per_L")
## # A tibble: 4 × 4
## pre_post_intervention variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 pre_Yellow nmol_per_L 0.894 0.157
## 2 post_Yellow nmol_per_L 0.861 0.0596
## 3 pre_Red nmol_per_L 0.962 0.815
## 4 post_Red nmol_per_L 0.867 0.0593
P > 0.05 for all 4 groups, suggesting data here is normal.
compare_means(nmol_per_L ~ pre_post, meta_table_lyc_long_rm_outlier, method = "t.test", paired = TRUE, group.by = "intervention")
## # A tibble: 2 × 9
## intervention .y. group2 group1 p p.adj p.format p.signif method
## <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
## 1 Red nmol_per_L post pre 0.00104 0.0021 0.001 ** T-test
## 2 Yellow nmol_per_L post pre 0.0934 0.093 0.093 ns T-test
(lyc_bp_stats_noOutlier <- meta_table_lyc_long_rm_outlier %>%
filter(intervention != "Baseline") %>%
ggpaired(x = "pre_post", y = "nmol_per_L", fill = "intervention",
line.color = "gray",
line.size = 1,
facet.by = "intervention",
short.panel.labs = FALSE, panel.labs = list(intervention = c("", "", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", linewidth = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "nmol/L plasma",
title = "Concentration of Lycopene",
subtitle = "") +
stat_compare_means(method = "t.test", paired = TRUE) )
Total lycopene levels are significantly increasing only after post-Red intervention
# calculate how many NAs there are per cytokine in whole data set
containsNAs_cytokines <- meta_table_edited %>%
filter(intervention != "Baseline") %>%
pivot_longer(cols = if_ng:il_4,
names_to = "cytokine",
values_to = "cyto_conc_pg_ml") %>%
group_by(cytokine, patient_id) %>%
count(is.na(cyto_conc_pg_ml)) %>%
filter(`is.na(cyto_conc_pg_ml)` == TRUE)
kable(containsNAs_cytokines)
| cytokine | patient_id | is.na(cyto_conc_pg_ml) | n |
|---|---|---|---|
| il_10 | 6101 | TRUE | 1 |
| il_10 | 6110 | TRUE | 2 |
| il_13 | 6103 | TRUE | 2 |
| il_13 | 6106 | TRUE | 3 |
| il_13 | 6109 | TRUE | 3 |
| il_13 | 6111 | TRUE | 3 |
| il_2 | 6103 | TRUE | 3 |
| il_2 | 6108 | TRUE | 1 |
| il_2 | 6109 | TRUE | 2 |
| il_2 | 6110 | TRUE | 4 |
| il_2 | 6112 | TRUE | 4 |
| il_4 | 6101 | TRUE | 4 |
| il_4 | 6111 | TRUE | 4 |
| il_4 | 6112 | TRUE | 4 |
| il_4 | 6113 | TRUE | 4 |
# calculate how many NAs there are per cytokine in whole data set
containsNAs_cells <- meta_table_edited %>%
filter(intervention != "Baseline") %>%
pivot_longer(cols = starts_with("x"),
names_to = "cell_type",
values_to = "cell_value") %>%
group_by(cell_type, pre_post_intervention) %>%
count(is.na(cell_value)) %>%
filter(`is.na(cell_value)` == TRUE)
kable(containsNAs_cells)
| cell_type | pre_post_intervention | is.na(cell_value) | n |
|---|
# impute any missing values by replacing them with 1/2 of the lowest concentration value of a cytokine (i.e. in a row).
imputed_meta_table <- meta_table_edited %>%
filter(intervention != "Baseline")
imputed_meta_table[,11:25] <- lapply(imputed_meta_table[,11:25],
function(x) ifelse(is.na(x),
min(x, na.rm = TRUE)/2, x))
dim(imputed_meta_table)
## [1] 48 93
# convert cytokines from wide to long
cytokines_long <- imputed_meta_table[-c(26:82)] %>%
pivot_longer(cols = if_ng:il_4,
names_to = "cytokine",
values_to = "cyto_conc_pg_ml")
Let’s remove subject data points for subject 6112’s yellow intervention first (their total lycopene levels increased post-yellow intervention)
cytokines_noOutlier_long <- anti_join(cytokines_long,
subj6112_yellow)
After testing several mixed linear models, the best model turned out to have pre_post as a fixed variable and patient_id as random effect. Carotenoids added no statistically significant effect to the model, suggesting that carotenoids may not contribute to cytokine levels. We decided to move on with paired mean comparisons for each group. May revisit mixed linear modeling to see if significant metabolites from metabolomics analyses affect immune outcomes (i.e. cytokine concentrations and immune cell populations).
cytokines_noOutlier_long %>%
group_by(cytokine) %>%
shapiro_test(cyto_conc_pg_ml)
## # A tibble: 15 × 4
## cytokine variable statistic p
## <chr> <chr> <dbl> <dbl>
## 1 gm_csf cyto_conc_pg_ml 0.613 8.76e-10
## 2 if_ng cyto_conc_pg_ml 0.596 4.99e-10
## 3 il_10 cyto_conc_pg_ml 0.736 9.45e- 8
## 4 il_12p40 cyto_conc_pg_ml 0.870 1.08e- 4
## 5 il_12p70 cyto_conc_pg_ml 0.416 2.74e-12
## 6 il_13 cyto_conc_pg_ml 0.691 1.48e- 8
## 7 il_1b cyto_conc_pg_ml 0.761 2.90e- 7
## 8 il_1ra cyto_conc_pg_ml 0.699 2.03e- 8
## 9 il_2 cyto_conc_pg_ml 0.552 1.23e-10
## 10 il_4 cyto_conc_pg_ml 0.641 2.29e- 9
## 11 il_5 cyto_conc_pg_ml 0.720 4.82e- 8
## 12 il_6 cyto_conc_pg_ml 0.652 3.42e- 9
## 13 il_8 cyto_conc_pg_ml 0.967 2.09e- 1
## 14 mcp_1 cyto_conc_pg_ml 0.931 9.13e- 3
## 15 tn_fa cyto_conc_pg_ml 0.676 8.57e- 9
Data is not normally distributed, so will use nonparametric paired tests for the mean comparisons.
First, let’s test for sequence effects. We don’t expect to have this problem since the participants did a 4-week washout period (no tomato/lycopene or soy isoflavone-containing foods). This should have been sufficient enough for the intervention to not have lingering effects on cytokine levels.
(wilcox_cyto_seqeffects <- cytokines_long %>%
filter(pre_post == "post") %>%
group_by(cytokine) %>%
wilcox_test(cyto_conc_pg_ml ~ sequence, paired = TRUE, p.adjust.method = "BH", detailed = TRUE))
## # A tibble: 15 × 13
## cytokine estimate .y. group1 group2 n1 n2 statistic p conf.low
## * <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 gm_csf 11.0 cyto_c… R_Y Y_R 12 12 52 0.339 -12.1
## 2 if_ng 0.638 cyto_c… R_Y Y_R 12 12 59 0.129 -0.815
## 3 il_10 -1.28 cyto_c… R_Y Y_R 12 12 15 0.12 -4.49
## 4 il_12p40 -1.92 cyto_c… R_Y Y_R 12 12 37 0.91 -33.9
## 5 il_12p70 0.19 cyto_c… R_Y Y_R 12 12 43 0.791 -1.66
## 6 il_13 18.3 cyto_c… R_Y Y_R 12 12 50 0.142 -4.31
## 7 il_1b 3.09 cyto_c… R_Y Y_R 12 12 49 0.47 -16.6
## 8 il_1ra 0.728 cyto_c… R_Y Y_R 12 12 48 0.519 -1.07
## 9 il_2 0.682 cyto_c… R_Y Y_R 12 12 40.5 0.533 -0.910
## 10 il_4 0.210 cyto_c… R_Y Y_R 12 12 33 0.61 -2.15
## 11 il_5 -1.19 cyto_c… R_Y Y_R 12 12 31 0.569 -4.70
## 12 il_6 0.832 cyto_c… R_Y Y_R 12 12 56.5 0.182 -0.435
## 13 il_8 -0.538 cyto_c… R_Y Y_R 12 12 29 0.47 -1.90
## 14 mcp_1 -32.5 cyto_c… R_Y Y_R 12 12 17 0.0923 -107.
## 15 tn_fa 4.97 cyto_c… R_Y Y_R 12 12 42 0.85 -8.93
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
For every cytokine, sequence does not significantly effect the outcome. Therefore we can continue to assume there are no sequence effects.
Non-parametric alternative to repeated measures ANOVA
#not working for me
friedman_test <- cytokines_long %>%
select(patient_id, pre_post_intervention, cytokine, cyto_conc_pg_ml) %>%
friedman_test(cyto_conc_pg_ml ~ pre_post_intervention | patient_id)
Comparing pre- to post-yellow (low carotenoid) intervention
(wilcox_cyto_yellow <- cytokines_noOutlier_long %>%
filter(intervention == "Yellow") %>%
group_by(cytokine) %>%
wilcox_test(cyto_conc_pg_ml ~ pre_post, paired = TRUE, p.adjust.method = "BH", detailed = TRUE))
## # A tibble: 15 × 13
## cytokine estimate .y. group1 group2 n1 n2 statistic p conf.low
## * <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 gm_csf 1.02 cyto_… pre post 11 11 42 0.465 -2.41
## 2 if_ng 0.0375 cyto_… pre post 11 11 41 0.52 -0.075
## 3 il_10 -0.0325 cyto_… pre post 11 11 32 0.966 -0.310
## 4 il_12p40 1.18 cyto_… pre post 11 11 29 0.919 -9.23
## 5 il_12p70 -0.110 cyto_… pre post 11 11 24 0.76 -0.340
## 6 il_13 -0.0000123 cyto_… pre post 11 11 22.5 1 -1.96
## 7 il_1b 3.73 cyto_… pre post 11 11 34 0.193 -0.965
## 8 il_1ra 0.348 cyto_… pre post 11 11 49 0.175 -0.340
## 9 il_2 0.285 cyto_… pre post 11 11 34 0.193 -0.0950
## 10 il_4 0.0535 cyto_… pre post 11 11 28 0.183 -0.0400
## 11 il_5 -0.538 cyto_… pre post 11 11 25 0.52 -2.01
## 12 il_6 0.868 cyto_… pre post 11 11 46 0.278 -0.805
## 13 il_8 -0.178 cyto_… pre post 11 11 27 0.638 -0.94
## 14 mcp_1 -1.72 cyto_… pre post 11 11 27 0.638 -16.2
## 15 tn_fa -0.373 cyto_… pre post 11 11 32 0.966 -2.93
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
# extract statistically significant cytokines
(sig_cytokines_yellow <- wilcox_cyto_yellow %>%
filter(p < 0.05))
## # A tibble: 0 × 13
## # ℹ 13 variables: cytokine <chr>, estimate <dbl>, .y. <chr>, group1 <chr>,
## # group2 <chr>, n1 <int>, n2 <int>, statistic <dbl>, p <dbl>, conf.low <dbl>,
## # conf.high <dbl>, method <chr>, alternative <chr>
Comparing pre- to post-red (tomato-soy) intervention
(wilcox_cyto_red <- cytokines_noOutlier_long %>%
filter(intervention == "Red") %>%
group_by(cytokine) %>%
wilcox_test(cyto_conc_pg_ml ~ pre_post, paired = TRUE, p.adjust.method = "BH", detailed = TRUE))
## # A tibble: 15 × 13
## cytokine estimate .y. group1 group2 n1 n2 statistic p conf.low
## * <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 gm_csf 2.20 cyto_c… pre post 12 12 65 0.0425 0.0350
## 2 if_ng 0.145 cyto_c… pre post 12 12 49 0.168 -0.0750
## 3 il_10 0.148 cyto_c… pre post 12 12 51 0.38 -0.12
## 4 il_12p40 10.0 cyto_c… pre post 12 12 61 0.0923 -2.34
## 5 il_12p70 0.672 cyto_c… pre post 12 12 66 0.0376 0.0401
## 6 il_13 1.67 cyto_c… pre post 12 12 47 0.23 -1.28
## 7 il_1b 0.986 cyto_c… pre post 12 12 48 0.197 -0.535
## 8 il_1ra 0.115 cyto_c… pre post 12 12 44 0.733 -0.405
## 9 il_2 0.170 cyto_c… pre post 12 12 31 0.343 -0.120
## 10 il_4 0.0922 cyto_c… pre post 12 12 26 0.294 -0.100
## 11 il_5 0.958 cyto_c… pre post 12 12 65 0.0425 0.0300
## 12 il_6 0.375 cyto_c… pre post 12 12 57 0.176 -0.165
## 13 il_8 -0.0765 cyto_c… pre post 12 12 31.5 0.583 -0.450
## 14 mcp_1 -0.875 cyto_c… pre post 12 12 38 0.97 -12.5
## 15 tn_fa 2.72 cyto_c… pre post 12 12 64 0.0522 -0.25
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
# extract statistically significant cytokines
(sig_cytokines_red <- wilcox_cyto_red %>%
filter(wilcox_cyto_red$p < 0.05))
## # A tibble: 3 × 13
## cytokine estimate .y. group1 group2 n1 n2 statistic p conf.low
## <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 gm_csf 2.20 cyto_co… pre post 12 12 65 0.0425 0.0350
## 2 il_12p70 0.672 cyto_co… pre post 12 12 66 0.0376 0.0401
## 3 il_5 0.958 cyto_co… pre post 12 12 65 0.0425 0.0300
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
cytokines_long %>%
filter(cytokine == "gm_csf") %>%
ggpaired(x = "pre_post", y = "cyto_conc_pg_ml", fill = "intervention", facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", "", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "pg/mL plasma",
title = "Concentration of GM-CSF",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
cytokines_long %>%
filter(cytokine == "gm_csf") %>%
filter(patient_id != 6102) %>% #6102 has really high levels so let's see what the data looks like without them
ggpaired(x = "pre_post", y = "cyto_conc_pg_ml", fill = "intervention", facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", "", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "pg/mL plasma",
title = "Concentration of GM-CSF",
subtitle = "")
meta_table %>%
filter(intervention != "Baseline") %>%
ggplot(aes(x = pre_post, y = gm_csf, color = intervention)) +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_classic() +
labs(x = "",
y = "GM-CSF conc (pg/mL)",
title = "GM-CSF levels in each patient pre- and post- red and yellow intervention")
cytokines_long %>%
filter(cytokine == "il_12p70") %>%
ggpaired(x = "pre_post", y = "cyto_conc_pg_ml", fill = "intervention", facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", "", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "pg/mL plasma",
title = "Concentration of IL-12p70",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
cytokines_long %>%
filter(cytokine == "il_12p70") %>%
filter(patient_id != 6102) %>% # take off subject 6102 since their levels are so high
ggpaired(x = "pre_post", y = "cyto_conc_pg_ml", fill = "intervention", facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", "", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "pg/mL plasma",
title = "Concentration of IL-12p70",
subtitle = "")
meta_table %>%
filter(intervention != "Baseline") %>%
ggplot(aes(x = pre_post, y = il_12p70, color = intervention)) +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_classic() +
labs(x = "",
y = "IL-12p70 conc (pg/mL)",
title = "IL-12p70 levels in each patient pre- and post- red and yellow intervention")
cytokines_long %>%
filter(cytokine == "il_5") %>%
filter(intervention != "Baseline") %>%
ggpaired(x = "pre_post", y = "cyto_conc_pg_ml", fill = "intervention", facet.by = "intervention", short.panel.labs = FALSE, panel.labs = list(intervention = c("", "", ""))) +
scale_fill_manual(values = c("Red" = "tomato1",
"Yellow" = "yellow1"),
labels = c("Control", "Tomato-Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "pg/mL plasma",
title = "Concentration of IL-5",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
meta_table %>%
filter(intervention != "Baseline") %>%
ggplot(aes(x = pre_post, y = il_5, color = intervention)) +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_classic() +
labs(x = "",
y = "IL-5 conc (pg/mL)",
title = "IL-5 levels in each patient pre- and post- red and yellow intervention")
Let’s look at the average fold change for significantly changing cells in yellow. We’ll look at the avg fold change comparing pre to post yellow, pre to post Red, and post-intervention comparison.
red_sigcytokines_subset <- cytokines_noOutlier_long %>%
filter(cytokine %in% sig_cytokines_red$cytokine)%>%
select(patient_id, pre_post_intervention, cytokine, cyto_conc_pg_ml) %>%
group_by(cytokine) %>%
pivot_wider(names_from = pre_post_intervention,
values_from = cyto_conc_pg_ml) %>%
mutate(yellow_FC = post_Yellow/pre_Yellow,
red_FC = post_Red/pre_Red,
post_intervention_FC = post_Red/post_Yellow)
red_sigcytokines_subset %>%
summarize(mean_yellow_FC = mean(yellow_FC, na.rm = TRUE),
mean_red_FC = mean(red_FC),
mean_intervention_FC = mean(post_intervention_FC, na.rm = TRUE))
## # A tibble: 3 × 4
## cytokine mean_yellow_FC mean_red_FC mean_intervention_FC
## <chr> <dbl> <dbl> <dbl>
## 1 gm_csf 0.954 0.830 0.916
## 2 il_12p70 1.02 0.735 0.916
## 3 il_5 1.32 0.828 0.817
(wilcox_cyto_intervention <- cytokines_noOutlier_long %>%
filter(patient_id != 6112) %>% # remove subject 6112 so the red vs yellow comparison is even
filter(pre_post == "post") %>%
group_by(cytokine) %>%
wilcox_test(cyto_conc_pg_ml ~ pre_post_intervention, paired = TRUE, p.adjust.method = "BH", detailed = TRUE))
## # A tibble: 15 × 13
## cytokine estimate .y. group1 group2 n1 n2 statistic p conf.low
## * <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 gm_csf -2.37 cyto_co… post_… post_… 11 11 30 0.831 -38.4
## 2 if_ng -0.122 cyto_co… post_… post_… 11 11 30 0.831 -1.6
## 3 il_10 0.193 cyto_co… post_… post_… 11 11 35 0.898 -2.51
## 4 il_12p40 2.99 cyto_co… post_… post_… 11 11 39 0.638 -35.8
## 5 il_12p70 0.46 cyto_co… post_… post_… 11 11 40 0.577 -6
## 6 il_13 -0.715 cyto_co… post_… post_… 11 11 27 1 -27.8
## 7 il_1b 0.885 cyto_co… post_… post_… 11 11 35 0.898 -11.6
## 8 il_1ra 0.64 cyto_co… post_… post_… 11 11 40 0.577 -1.78
## 9 il_2 -0.00998 cyto_co… post_… post_… 11 11 27 1 -1.96
## 10 il_4 -0.0513 cyto_co… post_… post_… 11 11 30 0.831 -1.84
## 11 il_5 1.34 cyto_co… post_… post_… 11 11 44 0.365 -3.29
## 12 il_6 0.788 cyto_co… post_… post_… 11 11 47 0.24 -0.98
## 13 il_8 -0.0150 cyto_co… post_… post_… 11 11 33 1 -1.08
## 14 mcp_1 -2.36 cyto_co… post_… post_… 11 11 32 0.966 -50.3
## 15 tn_fa 1.12 cyto_co… post_… post_… 11 11 36 0.831 -11.0
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
# extract statistically significant cytokines
(sig_cytokines_intervention <- wilcox_cyto_intervention %>%
filter(p < 0.05))
## # A tibble: 0 × 13
## # ℹ 13 variables: cytokine <chr>, estimate <dbl>, .y. <chr>, group1 <chr>,
## # group2 <chr>, n1 <int>, n2 <int>, statistic <dbl>, p <dbl>, conf.low <dbl>,
## # conf.high <dbl>, method <chr>, alternative <chr>
# convert immune cell data from wide to long
cells_long <- imputed_meta_table %>%
pivot_longer(cols = starts_with("x"),
names_to = "cell_type",
values_to = "cell_value")
Let’s remove subject data points for subject 6112’s yellow intervention first (their total lycopene levels increased post-yellow intervention)
cells_noOutlier_long <- anti_join(cells_long,
subj6112_yellow)
After testing several mixed linear models, the best model turned out to have pre_post as a fixed variable and patient_id as random effect. Carotenoids added no statistically significant effect to the model, suggesting that carotenoids are not contributing to cytokine levels. We decided to move on with paired mean comparisons for each group. May revisit mixed linear modeling to see if significant metabolites from metabolomics analyses affect immune outcomes (i.e. cytokine concentrations and immune cell populations).
cells_noOutlier_long %>%
group_by(cell_type) %>%
shapiro_test(cell_value)
## # A tibble: 39 × 4
## cell_type variable statistic p
## <chr> <chr> <dbl> <dbl>
## 1 x01_cd45_cd66b_lymph_dc_mono cell_value 0.718 0.0000000282
## 2 x02_cd45_cd66b_grans cell_value 0.693 0.0000000102
## 3 x03_cd3_cd45_cd3_t_cells cell_value 0.964 0.144
## 4 x04_tc_rgd_cd3_ab_t_cells cell_value 0.965 0.161
## 5 x05_cd4_cd8_cd8_t_cells cell_value 0.926 0.00490
## 6 x06_cd45ro_cd45ra_naive_cd8 cell_value 0.878 0.000135
## 7 x07_cd46ro_cd45ra_cm_cd8 cell_value 0.880 0.000151
## 8 x08_cd45ro_cd45ra_em_cd8 cell_value 0.892 0.000354
## 9 x09_cd45r0_cd45ra_te_cd8 cell_value 0.854 0.0000287
## 10 x10_cd38_hladr_activated_cd8 cell_value 0.960 0.106
## # ℹ 29 more rows
Data is not normally distributed for several cell types, so it will be more appropriate use nonparametric paired tests for the mean comparisons.
First, let’s test for sequence effects. We don’t expect to have this problem since the participants did a 4-week washout period (no tomato/lycopene or soy isoflavone-containing foods). This should have been sufficient enough for the intervention to not have lingering effects on cell populations.
wilcox_cells_seqeffects <- cells_long %>%
filter(pre_post == "post") %>%
group_by(cell_type) %>%
wilcox_test(cell_value ~ sequence, paired = TRUE, p.adjust.method = "BH", detailed = TRUE)
kable(wilcox_cells_seqeffects, format = "markdown", digits = 3)
| cell_type | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | 0.749 | cell_value | R_Y | Y_R | 12 | 12 | 46 | 0.622 | -1.927 | 3.530 | Wilcoxon | two.sided |
| x02_cd45_cd66b_grans | 0.104 | cell_value | R_Y | Y_R | 12 | 12 | 48 | 0.519 | -1.184 | 0.451 | Wilcoxon | two.sided |
| x03_cd3_cd45_cd3_t_cells | 0.547 | cell_value | R_Y | Y_R | 12 | 12 | 39 | 1.000 | -10.730 | 10.298 | Wilcoxon | two.sided |
| x04_tc_rgd_cd3_ab_t_cells | 0.852 | cell_value | R_Y | Y_R | 12 | 12 | 40 | 0.970 | -9.714 | 10.626 | Wilcoxon | two.sided |
| x05_cd4_cd8_cd8_t_cells | -3.533 | cell_value | R_Y | Y_R | 12 | 12 | 20 | 0.151 | -10.151 | 3.026 | Wilcoxon | two.sided |
| x06_cd45ro_cd45ra_naive_cd8 | -0.648 | cell_value | R_Y | Y_R | 12 | 12 | 33 | 0.677 | -4.287 | 2.635 | Wilcoxon | two.sided |
| x07_cd46ro_cd45ra_cm_cd8 | -0.092 | cell_value | R_Y | Y_R | 12 | 12 | 30 | 0.519 | -0.297 | 0.153 | Wilcoxon | two.sided |
| x08_cd45ro_cd45ra_em_cd8 | -1.191 | cell_value | R_Y | Y_R | 12 | 12 | 28 | 0.424 | -4.496 | 1.720 | Wilcoxon | two.sided |
| x09_cd45r0_cd45ra_te_cd8 | -0.080 | cell_value | R_Y | Y_R | 12 | 12 | 37 | 0.910 | -0.892 | 0.969 | Wilcoxon | two.sided |
| x10_cd38_hladr_activated_cd8 | -0.020 | cell_value | R_Y | Y_R | 12 | 12 | 30 | 0.519 | -0.072 | 0.035 | Wilcoxon | two.sided |
| x11_cd4_cd8_cd4_t_cells | 2.972 | cell_value | R_Y | Y_R | 12 | 12 | 49 | 0.470 | -6.548 | 13.344 | Wilcoxon | two.sided |
| x12_cd45ro_cd45ra_naive_cd4 | -0.139 | cell_value | R_Y | Y_R | 12 | 12 | 38 | 0.970 | -6.918 | 6.436 | Wilcoxon | two.sided |
| x13_cd45ro_cd45ra_cm_cd4 | -0.582 | cell_value | R_Y | Y_R | 12 | 12 | 32 | 0.622 | -2.699 | 8.300 | Wilcoxon | two.sided |
| x14_cd45ro_cd45ra | -0.957 | cell_value | R_Y | Y_R | 12 | 12 | 20 | 0.151 | -2.246 | 0.572 | Wilcoxon | two.sided |
| x15_cd45ro_cd45ra_te_cd4 | -0.223 | cell_value | R_Y | Y_R | 12 | 12 | 24 | 0.266 | -0.497 | 0.246 | Wilcoxon | two.sided |
| x16_cd38_hladr_activated_cd4 | -0.027 | cell_value | R_Y | Y_R | 12 | 12 | 32 | 0.622 | -0.234 | 0.178 | Wilcoxon | two.sided |
| x17_cd25_cd127_tregs | -0.137 | cell_value | R_Y | Y_R | 12 | 12 | 25 | 0.301 | -0.764 | 0.089 | Wilcoxon | two.sided |
| x18_ccr4_cd4_total_ccr4_treg | -0.214 | cell_value | R_Y | Y_R | 12 | 12 | 18 | 0.110 | -0.769 | 0.049 | Wilcoxon | two.sided |
| x19_cd45ra_cd45ro_ccr4_treg_naive | -0.019 | cell_value | R_Y | Y_R | 12 | 12 | 11 | 0.027 | -0.041 | -0.003 | Wilcoxon | two.sided |
| x20_hladr_total_ccr4_treg_activated | -0.048 | cell_value | R_Y | Y_R | 12 | 12 | 33 | 0.677 | -0.149 | 0.105 | Wilcoxon | two.sided |
| x21_cd45ra_cd45ro_ccr4_treg_memory | -0.130 | cell_value | R_Y | Y_R | 12 | 12 | 23 | 0.233 | -0.616 | 0.081 | Wilcoxon | two.sided |
| x22_cxcr3_ccr6_th1 | -0.975 | cell_value | R_Y | Y_R | 12 | 12 | 18 | 0.110 | -1.831 | 0.377 | Wilcoxon | two.sided |
| x23_cxcr3_ccr6_th2 | -2.984 | cell_value | R_Y | Y_R | 12 | 12 | 29 | 0.470 | -10.986 | 3.824 | Wilcoxon | two.sided |
| x24_cxcr3_ccr6_th17 | -2.054 | cell_value | R_Y | Y_R | 12 | 12 | 34 | 0.733 | -4.400 | 7.155 | Wilcoxon | two.sided |
| x25_cd19_cd3_b_cells | 1.007 | cell_value | R_Y | Y_R | 12 | 12 | 46 | 0.622 | -4.840 | 2.912 | Wilcoxon | two.sided |
| x26_cd27_ig_d_naive_b_cells | 0.299 | cell_value | R_Y | Y_R | 12 | 12 | 41 | 0.910 | -5.320 | 2.972 | Wilcoxon | two.sided |
| x27_cd27_ig_d_memory_b_cells | 0.335 | cell_value | R_Y | Y_R | 12 | 12 | 59 | 0.129 | -0.217 | 0.872 | Wilcoxon | two.sided |
| x28_cd27_ig_d_memory_resting_b_cells | 0.024 | cell_value | R_Y | Y_R | 12 | 12 | 42 | 0.850 | -0.165 | 0.546 | Wilcoxon | two.sided |
| x30_cd27_cd38_plasmablasts | 0.017 | cell_value | R_Y | Y_R | 12 | 12 | 51 | 0.380 | -0.009 | 0.066 | Wilcoxon | two.sided |
| x31_cd14_monocytes | 1.619 | cell_value | R_Y | Y_R | 12 | 12 | 45 | 0.677 | -6.241 | 12.224 | Wilcoxon | two.sided |
| x32_cd16_non_classical_mono | 0.171 | cell_value | R_Y | Y_R | 12 | 12 | 46 | 0.622 | -0.467 | 0.805 | Wilcoxon | two.sided |
| x33_cd16_classical_mono | 0.189 | cell_value | R_Y | Y_R | 12 | 12 | 44 | 0.733 | -7.760 | 11.174 | Wilcoxon | two.sided |
| x34_hladr_cd56 | 1.477 | cell_value | R_Y | Y_R | 12 | 12 | 42 | 0.850 | -5.279 | 9.547 | Wilcoxon | two.sided |
| x35_cd16_cd123_cd11c_p_dc | 0.020 | cell_value | R_Y | Y_R | 12 | 12 | 45 | 0.677 | -0.093 | 0.169 | Wilcoxon | two.sided |
| x36_cd16_cd123_cd11c_m_dc | 0.570 | cell_value | R_Y | Y_R | 12 | 12 | 40 | 0.970 | -7.313 | 9.070 | Wilcoxon | two.sided |
| x37_cd56_cd161_cd123_nk_cells | 0.975 | cell_value | R_Y | Y_R | 12 | 12 | 45 | 0.677 | -3.287 | 6.313 | Wilcoxon | two.sided |
| x38_cd16_nk_cells | -0.376 | cell_value | R_Y | Y_R | 12 | 12 | 36 | 0.850 | -3.960 | 2.063 | Wilcoxon | two.sided |
| x40_cd14_mdsc_mono | 0.522 | cell_value | R_Y | Y_R | 12 | 12 | 43 | 0.791 | -5.934 | 14.504 | Wilcoxon | two.sided |
| x41_cd66b_mdsc_grans | 0.013 | cell_value | R_Y | Y_R | 12 | 12 | 41 | 0.910 | -0.219 | 0.314 | Wilcoxon | two.sided |
# extract statistically significant cytokines
wilcox_cells_seqeffects %>%
filter(p < 0.05)
## # A tibble: 1 × 13
## cell_type estimate .y. group1 group2 n1 n2 statistic p conf.low
## <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 x19_cd45ra… -0.0193 cell… R_Y Y_R 12 12 11 0.0269 -0.0409
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
Seems to be sequence effects for cell type #19. We will keep this in mind for the rest of the analyses.
Comparing pre- to post-yellow (low carotenoid) intervention
wilcox_cells_yellow <- cells_noOutlier_long %>%
filter(intervention == "Yellow") %>%
group_by(cell_type) %>%
wilcox_test(cell_value ~ pre_post, paired = TRUE, p.adjust.method = "BH", detailed = TRUE)
kable(wilcox_cells_yellow, format = "markdown", digits = 3)
| cell_type | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | -0.020 | cell_value | pre | post | 12 | 12 | 39 | 1.000 | -2.423 | 1.338 | Wilcoxon | two.sided |
| x02_cd45_cd66b_grans | 0.157 | cell_value | pre | post | 12 | 12 | 57 | 0.176 | -0.044 | 0.609 | Wilcoxon | two.sided |
| x03_cd3_cd45_cd3_t_cells | 0.059 | cell_value | pre | post | 12 | 12 | 39 | 1.000 | -10.501 | 10.828 | Wilcoxon | two.sided |
| x04_tc_rgd_cd3_ab_t_cells | 0.458 | cell_value | pre | post | 12 | 12 | 42 | 0.850 | -10.329 | 11.068 | Wilcoxon | two.sided |
| x05_cd4_cd8_cd8_t_cells | -1.858 | cell_value | pre | post | 12 | 12 | 6 | 0.007 | -3.930 | -0.853 | Wilcoxon | two.sided |
| x06_cd45ro_cd45ra_naive_cd8 | -0.179 | cell_value | pre | post | 12 | 12 | 35 | 0.791 | -1.146 | 1.083 | Wilcoxon | two.sided |
| x07_cd46ro_cd45ra_cm_cd8 | 0.058 | cell_value | pre | post | 12 | 12 | 45 | 0.677 | -0.128 | 0.298 | Wilcoxon | two.sided |
| x08_cd45ro_cd45ra_em_cd8 | -1.252 | cell_value | pre | post | 12 | 12 | 0 | 0.000 | -1.989 | -0.695 | Wilcoxon | two.sided |
| x09_cd45r0_cd45ra_te_cd8 | -0.396 | cell_value | pre | post | 12 | 12 | 13 | 0.043 | -0.887 | -0.033 | Wilcoxon | two.sided |
| x10_cd38_hladr_activated_cd8 | -0.018 | cell_value | pre | post | 12 | 12 | 31 | 0.569 | -0.050 | 0.023 | Wilcoxon | two.sided |
| x11_cd4_cd8_cd4_t_cells | 3.011 | cell_value | pre | post | 12 | 12 | 49 | 0.470 | -5.959 | 11.409 | Wilcoxon | two.sided |
| x12_cd45ro_cd45ra_naive_cd4 | 1.816 | cell_value | pre | post | 12 | 12 | 52 | 0.339 | -2.224 | 6.323 | Wilcoxon | two.sided |
| x13_cd45ro_cd45ra_cm_cd4 | 1.004 | cell_value | pre | post | 12 | 12 | 48 | 0.519 | -1.876 | 3.551 | Wilcoxon | two.sided |
| x14_cd45ro_cd45ra | -0.444 | cell_value | pre | post | 12 | 12 | 17 | 0.092 | -0.974 | 0.182 | Wilcoxon | two.sided |
| x15_cd45ro_cd45ra_te_cd4 | -0.072 | cell_value | pre | post | 12 | 12 | 0 | 0.000 | -0.157 | -0.022 | Wilcoxon | two.sided |
| x16_cd38_hladr_activated_cd4 | -0.039 | cell_value | pre | post | 12 | 12 | 26 | 0.339 | -0.110 | 0.045 | Wilcoxon | two.sided |
| x17_cd25_cd127_tregs | -0.039 | cell_value | pre | post | 12 | 12 | 34 | 0.733 | -0.193 | 0.183 | Wilcoxon | two.sided |
| x18_ccr4_cd4_total_ccr4_treg | -0.055 | cell_value | pre | post | 12 | 12 | 29 | 0.470 | -0.184 | 0.125 | Wilcoxon | two.sided |
| x19_cd45ra_cd45ro_ccr4_treg_naive | -0.005 | cell_value | pre | post | 12 | 12 | 24 | 0.266 | -0.017 | 0.005 | Wilcoxon | two.sided |
| x20_hladr_total_ccr4_treg_activated | -0.039 | cell_value | pre | post | 12 | 12 | 25 | 0.301 | -0.107 | 0.073 | Wilcoxon | two.sided |
| x21_cd45ra_cd45ro_ccr4_treg_memory | -0.043 | cell_value | pre | post | 12 | 12 | 31 | 0.569 | -0.158 | 0.123 | Wilcoxon | two.sided |
| x22_cxcr3_ccr6_th1 | 0.308 | cell_value | pre | post | 12 | 12 | 49 | 0.470 | -0.481 | 1.096 | Wilcoxon | two.sided |
| x23_cxcr3_ccr6_th2 | 2.739 | cell_value | pre | post | 12 | 12 | 53 | 0.301 | -4.335 | 10.332 | Wilcoxon | two.sided |
| x24_cxcr3_ccr6_th17 | -0.330 | cell_value | pre | post | 12 | 12 | 28 | 0.424 | -4.462 | 0.484 | Wilcoxon | two.sided |
| x25_cd19_cd3_b_cells | -2.358 | cell_value | pre | post | 12 | 12 | 15 | 0.064 | -4.879 | 0.343 | Wilcoxon | two.sided |
| x26_cd27_ig_d_naive_b_cells | -2.076 | cell_value | pre | post | 12 | 12 | 11 | 0.027 | -3.807 | -0.062 | Wilcoxon | two.sided |
| x27_cd27_ig_d_memory_b_cells | -0.156 | cell_value | pre | post | 12 | 12 | 24 | 0.266 | -0.472 | 0.469 | Wilcoxon | two.sided |
| x28_cd27_ig_d_memory_resting_b_cells | 0.045 | cell_value | pre | post | 12 | 12 | 48 | 0.519 | -0.073 | 0.225 | Wilcoxon | two.sided |
| x30_cd27_cd38_plasmablasts | -0.001 | cell_value | pre | post | 12 | 12 | 36 | 0.850 | -0.038 | 0.015 | Wilcoxon | two.sided |
| x31_cd14_monocytes | 1.384 | cell_value | pre | post | 12 | 12 | 46 | 0.622 | -4.459 | 8.364 | Wilcoxon | two.sided |
| x32_cd16_non_classical_mono | 0.496 | cell_value | pre | post | 12 | 12 | 54 | 0.266 | -0.374 | 1.772 | Wilcoxon | two.sided |
| x33_cd16_classical_mono | 0.035 | cell_value | pre | post | 12 | 12 | 39 | 1.000 | -3.616 | 7.042 | Wilcoxon | two.sided |
| x34_hladr_cd56 | 0.578 | cell_value | pre | post | 12 | 12 | 40 | 0.970 | -4.647 | 6.576 | Wilcoxon | two.sided |
| x35_cd16_cd123_cd11c_p_dc | -0.010 | cell_value | pre | post | 12 | 12 | 35 | 0.791 | -0.094 | 0.073 | Wilcoxon | two.sided |
| x36_cd16_cd123_cd11c_m_dc | 0.346 | cell_value | pre | post | 12 | 12 | 42 | 0.850 | -3.829 | 6.239 | Wilcoxon | two.sided |
| x37_cd56_cd161_cd123_nk_cells | -0.565 | cell_value | pre | post | 12 | 12 | 5 | 0.005 | -1.385 | -0.172 | Wilcoxon | two.sided |
| x38_cd16_nk_cells | 0.169 | cell_value | pre | post | 12 | 12 | 41 | 0.910 | -1.645 | 4.287 | Wilcoxon | two.sided |
| x40_cd14_mdsc_mono | 0.051 | cell_value | pre | post | 12 | 12 | 53 | 0.301 | -0.049 | 0.219 | Wilcoxon | two.sided |
| x41_cd66b_mdsc_grans | 0.022 | cell_value | pre | post | 12 | 12 | 54 | 0.266 | -0.017 | 0.119 | Wilcoxon | two.sided |
# extract statistically significant cytokines
(sig_cells_yellow <- wilcox_cells_yellow %>%
filter(p < 0.05))
## # A tibble: 6 × 13
## cell_type estimate .y. group1 group2 n1 n2 statistic p conf.low
## <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 x05_cd4_c… -1.86 cell… pre post 12 12 6 6.84e-3 -3.93
## 2 x08_cd45r… -1.25 cell… pre post 12 12 0 4.88e-4 -1.99
## 3 x09_cd45r… -0.396 cell… pre post 12 12 13 4.25e-2 -0.887
## 4 x15_cd45r… -0.0723 cell… pre post 12 12 0 4.88e-4 -0.157
## 5 x26_cd27_… -2.08 cell… pre post 12 12 11 2.69e-2 -3.81
## 6 x37_cd56_… -0.565 cell… pre post 12 12 5 4.88e-3 -1.38
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
(sigcells_yellow_bp <- cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_yellow$cell_type) %>%
filter(intervention == "Yellow") %>%
ggpaired(x = "pre_post", y = "cell_value", fill = "intervention", facet.by = c("cell_type")) +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Control"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Cell percentage",
title = "Cell populations significantly different between pre- and post-Yellow",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH"))
Log2 the scale so change is more visible
cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_yellow$cell_type) %>%
filter(intervention == "Yellow") %>%
mutate(log2_cellvalue = log2(cell_value)) %>%
ggpaired(x = "pre_post", y = "log2_cellvalue", fill = "intervention", facet.by = c("cell_type")) +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Control"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Log2 cell %",
title = "Cell populations significantly different between pre- and post-Yellow",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
Let’s see what things look like in red intervention
cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_yellow$cell_type) %>%
mutate(log2_cellvalue = log2(cell_value)) %>%
filter(intervention == "Red") %>%
ggpaired(x = "pre_post", y = "log2_cellvalue", fill = "intervention", facet.by = c("cell_type")) +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Tomato Soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Log2 cell %",
title = "Cell populations significantly different between pre- and post-Yellow",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
Cell #26 (Naive B-cell) is significantly increased after both interventions, although post-red has a stronger significance.
To get a sense of cell population percentages on the individual basis, we’ll look at line graphs.
I’m going to make a function for this for this first.
# function for line plots
line_plots_cells_fx <- function(immune_cell) {
cells_long %>% # even though this subject 6112 wasnt included for yellow stats, I'll still include them for the plots
filter(cell_type == immune_cell) %>%
ggplot(aes(x = pre_post, y = cell_value, color = intervention)) +
geom_line(aes(group = intervention)) +
scale_color_manual(values = c("Yellow" = "gold",
"Red" = "tomato1")) +
facet_wrap(vars(patient_id), scales = "free_y") +
theme_classic() +
labs(x = "",
y = "Cell population (%)",
title = paste(immune_cell, "in each patient pre- and post- red and yellow intervention"))
}
# create list for plots
lineplots_cells <- list()
# run this function on all cell types
for (i in cells_noOutlier_long$cell_type) {
lineplots_cells[[i]] = line_plots_cells_fx(i)
}
# signif immune cells pre vs post yellow
sig_cells_yellow$cell_type
## [1] "x05_cd4_cd8_cd8_t_cells" "x08_cd45ro_cd45ra_em_cd8"
## [3] "x09_cd45r0_cd45ra_te_cd8" "x15_cd45ro_cd45ra_te_cd4"
## [5] "x26_cd27_ig_d_naive_b_cells" "x37_cd56_cd161_cd123_nk_cells"
lineplots_cells[sig_cells_yellow$cell_type[1]]
## $x05_cd4_cd8_cd8_t_cells
lineplots_cells[sig_cells_yellow$cell_type[2]]
## $x08_cd45ro_cd45ra_em_cd8
lineplots_cells[sig_cells_yellow$cell_type[3]]
## $x09_cd45r0_cd45ra_te_cd8
lineplots_cells[sig_cells_yellow$cell_type[4]]
## $x15_cd45ro_cd45ra_te_cd4
Let’s look at the average fold change for significantly changing cells in yellow. We’ll look at the avg fold change comparing pre to post yellow, pre to post Red, and post-intervention comparison.
yellow_sigcells_subset <- cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_yellow$cell_type)%>%
select(patient_id, pre_post_intervention, cell_type, cell_value) %>%
group_by(cell_type) %>%
pivot_wider(names_from = pre_post_intervention,
values_from = cell_value) %>%
mutate(yellow_FC = post_Yellow/pre_Yellow,
red_FC = post_Red/pre_Red,
post_intervention_FC = post_Red/post_Yellow)
yellow_sigcells_subset %>%
summarize(mean_yellow_FC = mean(yellow_FC, na.rm = TRUE),
mean_red_FC = mean(red_FC),
mean_intervention_FC = mean(post_intervention_FC, na.rm = TRUE))
## # A tibble: 6 × 4
## cell_type mean_yellow_FC mean_red_FC mean_intervention_FC
## <chr> <dbl> <dbl> <dbl>
## 1 x05_cd4_cd8_cd8_t_cells 1.33 1.31 1.06
## 2 x08_cd45ro_cd45ra_em_cd8 1.52 2.06 1.75
## 3 x09_cd45r0_cd45ra_te_cd8 1.36 1.23 0.972
## 4 x15_cd45ro_cd45ra_te_cd4 1.60 48.0 20.5
## 5 x26_cd27_ig_d_naive_b_cells 1.62 1.84 0.987
## 6 x37_cd56_cd161_cd123_nk_cells 1.52 1.20 3.62
# the fold change for cell type #15 was huge for red FC! subject 6108 had a very large increase post-red, so let's see what the fold change is like without that subject
yellow_sigcells_subset %>%
filter(patient_id != 6108) %>%
summarize(mean_yellow_FC = mean(yellow_FC, na.rm = TRUE),
mean_red_FC = mean(red_FC),
mean_intervention_FC = mean(post_intervention_FC, na.rm = TRUE))
## # A tibble: 6 × 4
## cell_type mean_yellow_FC mean_red_FC mean_intervention_FC
## <chr> <dbl> <dbl> <dbl>
## 1 x05_cd4_cd8_cd8_t_cells 1.32 1.10 0.933
## 2 x08_cd45ro_cd45ra_em_cd8 1.51 1.22 0.940
## 3 x09_cd45r0_cd45ra_te_cd8 1.35 1.28 1.01
## 4 x15_cd45ro_cd45ra_te_cd4 1.54 1.42 0.930
## 5 x26_cd27_ig_d_naive_b_cells 1.59 1.77 0.935
## 6 x37_cd56_cd161_cd123_nk_cells 1.44 1.30 3.93
Comparing pre- to post-yellow (low carotenoid) intervention
wilcox_cells_red <- cells_noOutlier_long %>%
filter(intervention == "Red") %>%
group_by(cell_type) %>%
wilcox_test(cell_value ~ pre_post, paired = TRUE, p.adjust.method = "BH", detailed = TRUE)
kable(wilcox_cells_red, format = "markdown", digits = 3)
| cell_type | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | 0.676 | cell_value | pre | post | 12 | 12 | 63 | 0.064 | -0.051 | 2.407 | Wilcoxon | two.sided |
| x02_cd45_cd66b_grans | -0.250 | cell_value | pre | post | 12 | 12 | 22 | 0.204 | -0.702 | 0.097 | Wilcoxon | two.sided |
| x03_cd3_cd45_cd3_t_cells | 0.392 | cell_value | pre | post | 12 | 12 | 40 | 0.970 | -12.266 | 14.584 | Wilcoxon | two.sided |
| x04_tc_rgd_cd3_ab_t_cells | 0.624 | cell_value | pre | post | 12 | 12 | 41 | 0.910 | -12.085 | 14.719 | Wilcoxon | two.sided |
| x05_cd4_cd8_cd8_t_cells | -1.134 | cell_value | pre | post | 12 | 12 | 17 | 0.092 | -3.764 | 0.352 | Wilcoxon | two.sided |
| x06_cd45ro_cd45ra_naive_cd8 | -0.206 | cell_value | pre | post | 12 | 12 | 33 | 0.677 | -1.184 | 0.832 | Wilcoxon | two.sided |
| x07_cd46ro_cd45ra_cm_cd8 | -0.047 | cell_value | pre | post | 12 | 12 | 26 | 0.339 | -0.176 | 0.054 | Wilcoxon | two.sided |
| x08_cd45ro_cd45ra_em_cd8 | -0.904 | cell_value | pre | post | 12 | 12 | 13 | 0.043 | -2.261 | -0.013 | Wilcoxon | two.sided |
| x09_cd45r0_cd45ra_te_cd8 | -0.130 | cell_value | pre | post | 12 | 12 | 31 | 0.569 | -0.562 | 0.268 | Wilcoxon | two.sided |
| x10_cd38_hladr_activated_cd8 | 0.003 | cell_value | pre | post | 12 | 12 | 44 | 0.733 | -0.032 | 0.039 | Wilcoxon | two.sided |
| x11_cd4_cd8_cd4_t_cells | 2.054 | cell_value | pre | post | 12 | 12 | 44 | 0.733 | -9.442 | 14.299 | Wilcoxon | two.sided |
| x12_cd45ro_cd45ra_naive_cd4 | 0.851 | cell_value | pre | post | 12 | 12 | 44 | 0.733 | -2.927 | 8.874 | Wilcoxon | two.sided |
| x13_cd45ro_cd45ra_cm_cd4 | 0.850 | cell_value | pre | post | 12 | 12 | 49 | 0.470 | -1.505 | 2.827 | Wilcoxon | two.sided |
| x14_cd45ro_cd45ra | -0.421 | cell_value | pre | post | 12 | 12 | 19 | 0.129 | -0.910 | 0.237 | Wilcoxon | two.sided |
| x15_cd45ro_cd45ra_te_cd4 | -0.052 | cell_value | pre | post | 12 | 12 | 18 | 0.110 | -0.178 | 0.008 | Wilcoxon | two.sided |
| x16_cd38_hladr_activated_cd4 | 0.028 | cell_value | pre | post | 12 | 12 | 52 | 0.339 | -0.025 | 0.121 | Wilcoxon | two.sided |
| x17_cd25_cd127_tregs | 0.025 | cell_value | pre | post | 12 | 12 | 50 | 0.424 | -0.100 | 0.219 | Wilcoxon | two.sided |
| x18_ccr4_cd4_total_ccr4_treg | 0.028 | cell_value | pre | post | 12 | 12 | 45 | 0.677 | -0.077 | 0.215 | Wilcoxon | two.sided |
| x19_cd45ra_cd45ro_ccr4_treg_naive | 0.002 | cell_value | pre | post | 12 | 12 | 45 | 0.677 | -0.007 | 0.014 | Wilcoxon | two.sided |
| x20_hladr_total_ccr4_treg_activated | 0.006 | cell_value | pre | post | 12 | 12 | 42 | 0.850 | -0.054 | 0.066 | Wilcoxon | two.sided |
| x21_cd45ra_cd45ro_ccr4_treg_memory | 0.037 | cell_value | pre | post | 12 | 12 | 48 | 0.519 | -0.064 | 0.157 | Wilcoxon | two.sided |
| x22_cxcr3_ccr6_th1 | -0.295 | cell_value | pre | post | 12 | 12 | 28 | 0.424 | -0.995 | 0.206 | Wilcoxon | two.sided |
| x23_cxcr3_ccr6_th2 | 2.253 | cell_value | pre | post | 12 | 12 | 51 | 0.380 | -3.040 | 15.172 | Wilcoxon | two.sided |
| x24_cxcr3_ccr6_th17 | 0.070 | cell_value | pre | post | 12 | 12 | 41 | 0.910 | -0.905 | 1.131 | Wilcoxon | two.sided |
| x25_cd19_cd3_b_cells | -3.122 | cell_value | pre | post | 12 | 12 | 5 | 0.005 | -4.960 | -0.788 | Wilcoxon | two.sided |
| x26_cd27_ig_d_naive_b_cells | -2.445 | cell_value | pre | post | 12 | 12 | 9 | 0.016 | -4.543 | -0.476 | Wilcoxon | two.sided |
| x27_cd27_ig_d_memory_b_cells | -0.174 | cell_value | pre | post | 12 | 12 | 24 | 0.266 | -0.511 | 0.198 | Wilcoxon | two.sided |
| x28_cd27_ig_d_memory_resting_b_cells | -0.061 | cell_value | pre | post | 12 | 12 | 16 | 0.077 | -0.130 | 0.008 | Wilcoxon | two.sided |
| x30_cd27_cd38_plasmablasts | 0.007 | cell_value | pre | post | 12 | 12 | 49 | 0.470 | -0.016 | 0.048 | Wilcoxon | two.sided |
| x31_cd14_monocytes | 1.280 | cell_value | pre | post | 12 | 12 | 43 | 0.791 | -7.902 | 12.256 | Wilcoxon | two.sided |
| x32_cd16_non_classical_mono | 0.578 | cell_value | pre | post | 12 | 12 | 49 | 0.470 | -0.695 | 1.786 | Wilcoxon | two.sided |
| x33_cd16_classical_mono | 0.725 | cell_value | pre | post | 12 | 12 | 39 | 1.000 | -7.526 | 9.846 | Wilcoxon | two.sided |
| x34_hladr_cd56 | 1.101 | cell_value | pre | post | 12 | 12 | 43 | 0.791 | -7.000 | 10.814 | Wilcoxon | two.sided |
| x35_cd16_cd123_cd11c_p_dc | 0.076 | cell_value | pre | post | 12 | 12 | 51 | 0.380 | -0.078 | 0.217 | Wilcoxon | two.sided |
| x36_cd16_cd123_cd11c_m_dc | 0.267 | cell_value | pre | post | 12 | 12 | 39 | 1.000 | -6.094 | 7.975 | Wilcoxon | two.sided |
| x37_cd56_cd161_cd123_nk_cells | 0.186 | cell_value | pre | post | 12 | 12 | 43 | 0.791 | -2.056 | 4.137 | Wilcoxon | two.sided |
| x38_cd16_nk_cells | 0.217 | cell_value | pre | post | 12 | 12 | 49 | 0.470 | -0.455 | 1.966 | Wilcoxon | two.sided |
| x40_cd14_mdsc_mono | 0.805 | cell_value | pre | post | 12 | 12 | 48 | 0.519 | -3.918 | 6.494 | Wilcoxon | two.sided |
| x41_cd66b_mdsc_grans | 0.005 | cell_value | pre | post | 12 | 12 | 43 | 0.791 | -0.232 | 0.086 | Wilcoxon | two.sided |
# extract statistically significant cytokines
(sig_cells_red <- wilcox_cells_red %>%
filter(p < 0.05))
## # A tibble: 3 × 13
## cell_type estimate .y. group1 group2 n1 n2 statistic p conf.low
## <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 x08_cd45r… -0.904 cell… pre post 12 12 13 0.0425 -2.26
## 2 x25_cd19_… -3.12 cell… pre post 12 12 5 0.00488 -4.96
## 3 x26_cd27_… -2.45 cell… pre post 12 12 9 0.0161 -4.54
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_red$cell_type) %>%
filter(intervention == "Red") %>%
ggpaired(x = "pre_post", y = "cell_value", fill = "intervention", facet.by = "cell_type", short.panel.labs = FALSE, panel.labs = list(cell_type = c("CD45RO+ CD45RA- (EM CD8)","CD19+CD3- (B cells)", "CD27-IgD+ (Naive B cells)"))) +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Tomato-soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Cell value",
title = "Cell populations significantly different between pre- and post-Tomato soy",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
Let’s see what this looks like pre- and post- control
cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_red$cell_type) %>%
filter(intervention == "Yellow") %>%
ggpaired(x = "pre_post", y = "cell_value", fill = "intervention", facet.by = "cell_type", short.panel.labs = FALSE, panel.labs = list(cell_type = c("CD45RO+ CD45RA- (EM CD8)","CD19+CD3- (B cells)", "CD27-IgD+ (Naive B cells)"))) +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Control"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.15) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Cell value",
title = "Cell populations significantly different between pre- and post-Tomato soy",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
EM CD8 and Naive B- cells are significantly increased post control as well. Stronger effect in control.
Let’s look at the average fold change for significantly changing cells in red. We’ll look at the avg fold change comparing pre to post yellow, pre to post Red, and post-intervention comparison.
red_sigcells_subset <- cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_red$cell_type)%>%
select(patient_id, pre_post_intervention, cell_type, cell_value) %>%
group_by(cell_type) %>%
pivot_wider(names_from = pre_post_intervention,
values_from = cell_value) %>%
mutate(yellow_FC = post_Yellow/pre_Yellow,
red_FC = post_Red/pre_Red,
post_intervention_FC = post_Red/post_Yellow)
red_sigcells_subset %>%
summarize(mean_yellow_FC = mean(yellow_FC, na.rm = TRUE),
mean_red_FC = mean(red_FC),
mean_intervention_FC = mean(post_intervention_FC, na.rm = TRUE))
## # A tibble: 3 × 4
## cell_type mean_yellow_FC mean_red_FC mean_intervention_FC
## <chr> <dbl> <dbl> <dbl>
## 1 x08_cd45ro_cd45ra_em_cd8 1.52 2.06 1.75
## 2 x25_cd19_cd3_b_cells 1.49 1.64 0.972
## 3 x26_cd27_ig_d_naive_b_cells 1.62 1.84 0.987
Comparing post-tomato soy to post-control intervention
wilcox_cells_intervention <- cells_noOutlier_long %>%
filter(pre_post == "post",
patient_id != "6112") %>%
group_by(cell_type) %>%
wilcox_test(cell_value ~ intervention, paired = TRUE, p.adjust.method = "BH", detailed = TRUE)
kable(wilcox_cells_intervention, format = "markdown", digits = 3)
| cell_type | estimate | .y. | group1 | group2 | n1 | n2 | statistic | p | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| x01_cd45_cd66b_lymph_dc_mono | 0.939 | cell_value | Yellow | Red | 11 | 11 | 40 | 0.577 | -2.155 | 4.409 | Wilcoxon | two.sided |
| x02_cd45_cd66b_grans | -0.332 | cell_value | Yellow | Red | 11 | 11 | 19 | 0.240 | -1.196 | 0.487 | Wilcoxon | two.sided |
| x03_cd3_cd45_cd3_t_cells | -6.112 | cell_value | Yellow | Red | 11 | 11 | 25 | 0.520 | -16.628 | 9.710 | Wilcoxon | two.sided |
| x04_tc_rgd_cd3_ab_t_cells | -5.992 | cell_value | Yellow | Red | 11 | 11 | 24 | 0.465 | -16.693 | 9.095 | Wilcoxon | two.sided |
| x05_cd4_cd8_cd8_t_cells | -1.102 | cell_value | Yellow | Red | 11 | 11 | 29 | 0.765 | -4.900 | 4.591 | Wilcoxon | two.sided |
| x06_cd45ro_cd45ra_naive_cd8 | -0.099 | cell_value | Yellow | Red | 11 | 11 | 33 | 1.000 | -3.458 | 3.790 | Wilcoxon | two.sided |
| x07_cd46ro_cd45ra_cm_cd8 | -0.023 | cell_value | Yellow | Red | 11 | 11 | 30 | 0.831 | -0.273 | 0.239 | Wilcoxon | two.sided |
| x08_cd45ro_cd45ra_em_cd8 | 0.085 | cell_value | Yellow | Red | 11 | 11 | 34 | 0.966 | -4.807 | 3.087 | Wilcoxon | two.sided |
| x09_cd45r0_cd45ra_te_cd8 | 0.353 | cell_value | Yellow | Red | 11 | 11 | 41 | 0.520 | -0.679 | 1.444 | Wilcoxon | two.sided |
| x10_cd38_hladr_activated_cd8 | 0.023 | cell_value | Yellow | Red | 11 | 11 | 50 | 0.147 | -0.013 | 0.060 | Wilcoxon | two.sided |
| x11_cd4_cd8_cd4_t_cells | -3.659 | cell_value | Yellow | Red | 11 | 11 | 24 | 0.465 | -14.703 | 6.861 | Wilcoxon | two.sided |
| x12_cd45ro_cd45ra_naive_cd4 | -1.463 | cell_value | Yellow | Red | 11 | 11 | 27 | 0.638 | -9.292 | 4.505 | Wilcoxon | two.sided |
| x13_cd45ro_cd45ra_cm_cd4 | -0.406 | cell_value | Yellow | Red | 11 | 11 | 30 | 0.831 | -5.530 | 5.506 | Wilcoxon | two.sided |
| x14_cd45ro_cd45ra | 0.190 | cell_value | Yellow | Red | 11 | 11 | 35 | 0.898 | -1.565 | 1.691 | Wilcoxon | two.sided |
| x15_cd45ro_cd45ra_te_cd4 | -0.006 | cell_value | Yellow | Red | 11 | 11 | 31 | 0.898 | -6.261 | 0.342 | Wilcoxon | two.sided |
| x16_cd38_hladr_activated_cd4 | 0.045 | cell_value | Yellow | Red | 11 | 11 | 43 | 0.413 | -0.118 | 0.217 | Wilcoxon | two.sided |
| x17_cd25_cd127_tregs | -0.030 | cell_value | Yellow | Red | 11 | 11 | 31 | 0.898 | -0.225 | 0.779 | Wilcoxon | two.sided |
| x18_ccr4_cd4_total_ccr4_treg | -0.010 | cell_value | Yellow | Red | 11 | 11 | 32 | 0.966 | -0.227 | 0.774 | Wilcoxon | two.sided |
| x19_cd45ra_cd45ro_ccr4_treg_naive | 0.004 | cell_value | Yellow | Red | 11 | 11 | 39 | 0.638 | -0.014 | 0.032 | Wilcoxon | two.sided |
| x20_hladr_total_ccr4_treg_activated | 0.071 | cell_value | Yellow | Red | 11 | 11 | 46 | 0.278 | -0.060 | 0.212 | Wilcoxon | two.sided |
| x21_cd45ra_cd45ro_ccr4_treg_memory | -0.012 | cell_value | Yellow | Red | 11 | 11 | 31 | 0.898 | -0.213 | 0.501 | Wilcoxon | two.sided |
| x22_cxcr3_ccr6_th1 | -0.319 | cell_value | Yellow | Red | 11 | 11 | 29 | 0.765 | -1.790 | 1.121 | Wilcoxon | two.sided |
| x23_cxcr3_ccr6_th2 | -2.599 | cell_value | Yellow | Red | 11 | 11 | 27 | 0.638 | -10.608 | 7.891 | Wilcoxon | two.sided |
| x24_cxcr3_ccr6_th17 | 0.108 | cell_value | Yellow | Red | 11 | 11 | 34 | 0.966 | -3.746 | 4.798 | Wilcoxon | two.sided |
| x25_cd19_cd3_b_cells | -0.635 | cell_value | Yellow | Red | 11 | 11 | 28 | 0.700 | -5.464 | 6.970 | Wilcoxon | two.sided |
| x26_cd27_ig_d_naive_b_cells | -0.199 | cell_value | Yellow | Red | 11 | 11 | 32 | 0.966 | -5.914 | 5.040 | Wilcoxon | two.sided |
| x27_cd27_ig_d_memory_b_cells | 0.191 | cell_value | Yellow | Red | 11 | 11 | 41 | 0.520 | -0.497 | 0.844 | Wilcoxon | two.sided |
| x28_cd27_ig_d_memory_resting_b_cells | 0.043 | cell_value | Yellow | Red | 11 | 11 | 34 | 0.966 | -0.341 | 0.346 | Wilcoxon | two.sided |
| x30_cd27_cd38_plasmablasts | 0.018 | cell_value | Yellow | Red | 11 | 11 | 43 | 0.413 | -0.020 | 0.067 | Wilcoxon | two.sided |
| x31_cd14_monocytes | 5.496 | cell_value | Yellow | Red | 11 | 11 | 42 | 0.465 | -7.390 | 12.925 | Wilcoxon | two.sided |
| x32_cd16_non_classical_mono | -0.170 | cell_value | Yellow | Red | 11 | 11 | 26 | 0.577 | -0.805 | 0.372 | Wilcoxon | two.sided |
| x33_cd16_classical_mono | 4.795 | cell_value | Yellow | Red | 11 | 11 | 43 | 0.413 | -7.726 | 13.178 | Wilcoxon | two.sided |
| x34_hladr_cd56 | 4.297 | cell_value | Yellow | Red | 11 | 11 | 41 | 0.520 | -6.792 | 11.745 | Wilcoxon | two.sided |
| x35_cd16_cd123_cd11c_p_dc | 0.070 | cell_value | Yellow | Red | 11 | 11 | 46 | 0.278 | -0.067 | 0.157 | Wilcoxon | two.sided |
| x36_cd16_cd123_cd11c_m_dc | 3.045 | cell_value | Yellow | Red | 11 | 11 | 42 | 0.465 | -6.255 | 12.109 | Wilcoxon | two.sided |
| x37_cd56_cd161_cd123_nk_cells | -4.028 | cell_value | Yellow | Red | 11 | 11 | 14 | 0.102 | -9.791 | 0.180 | Wilcoxon | two.sided |
| x38_cd16_nk_cells | 2.788 | cell_value | Yellow | Red | 11 | 11 | 63 | 0.005 | 1.068 | 5.805 | Wilcoxon | two.sided |
| x40_cd14_mdsc_mono | -7.563 | cell_value | Yellow | Red | 11 | 11 | 0 | 0.001 | -20.535 | -2.532 | Wilcoxon | two.sided |
| x41_cd66b_mdsc_grans | -0.011 | cell_value | Yellow | Red | 11 | 11 | 29 | 0.765 | -0.345 | 0.219 | Wilcoxon | two.sided |
# extract statistically significant cytokines
(sig_cells_intervention <- wilcox_cells_intervention %>%
filter(p < 0.05))
## # A tibble: 2 × 13
## cell_type estimate .y. group1 group2 n1 n2 statistic p conf.low
## <chr> <dbl> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 x38_cd16_… 2.79 cell… Yellow Red 11 11 63 4.88e-3 1.07
## 2 x40_cd14_… -7.56 cell… Yellow Red 11 11 0 9.77e-4 -20.5
## # ℹ 3 more variables: conf.high <dbl>, method <chr>, alternative <chr>
cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_intervention$cell_type) %>%
filter(pre_post == "post") %>%
ggpaired(x = "pre_post_intervention", y = "cell_value", fill = "intervention", facet.by = "cell_type") +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Control", "Tomato-soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.25) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Cell value",
title = "Cell populations significantly different between post-Control and post-Tomato soy juices",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH", comparisons = list(c("post_Yellow", "post_Red")), label.y = 40)
Let’s look at the trends for these cell types during individual interventions
cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_intervention$cell_type) %>%
filter(intervention == "Red") %>%
ggpaired(x = "pre_post", y = "cell_value", fill = "intervention", facet.by = "cell_type", short.panel.labs = FALSE, panel.labs = list(cell_type = c("CD16- NK cells","CD14+ MDSC (Mono)"))) +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Tomato-soy"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.25) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Cell value",
title = "Cell populations significantly different between post-Control and post-Tomato soy juices",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_intervention$cell_type) %>%
filter(intervention == "Yellow") %>%
ggpaired(x = "pre_post", y = "cell_value", fill = "intervention", facet.by = "cell_type", short.panel.labs = FALSE, panel.labs = list(cell_type = c("CD16- NK cells","CD14+ MDSC (Mono)"))) +
scale_fill_manual(values = c("Yellow" = "yellow1",
"Red" = "tomato1"),
labels = c("Control"),
name = "Intervention") +
geom_line(aes(group = patient_id), colour = "gray", size = 0.25) +
theme_clean(base_size = 18, base_family = "sans") +
labs(x = "",
y = "Cell value",
title = "Cell populations significantly different between post-Control and post-Tomato soy juices",
subtitle = "") +
stat_compare_means(method = "wilcox.test", paired = TRUE, p.adjust.method = "BH")
To get a sense of cell population percentages on the individual basis, we’ll look at line graphs.
lineplots_cells[sig_cells_intervention$cell_type[1]]
## $x38_cd16_nk_cells
Looking at post control v post tomato soy
imputed_meta_table %>%
filter(pre_post == "post") %>%
ggplot(aes(x = pre_post_intervention, y = x38_cd16_nk_cells, color = patient_id)) +
geom_line(aes(group = patient_id)) +
theme_bw() +
labs(x = "Timepoint",
y = "Population level (%)",
title = "CD16- NK Cell populations in each patient between post- interventions")
lineplots_cells[sig_cells_intervention$cell_type[2]]
## $x40_cd14_mdsc_mono
Let’s look at post control vs post tomato soy
imputed_meta_table %>%
filter(pre_post == "post") %>%
ggplot(aes(x = pre_post_intervention, y = x40_cd14_mdsc_mono, color = patient_id)) +
geom_line(aes(group = patient_id)) +
theme_bw() +
labs(x = "Timepoint",
y = "Population level (%)",
title = "CD14+ monocytic MDSC populations in each patient between post- interventions")
Let’s look at the average fold change for significantly changing cells in red. We’ll look at the avg fold change comparing pre to post yellow, pre to post Red, and post-intervention comparison.
intervention_sigcells_subset <- cells_noOutlier_long %>%
filter(cell_type %in% sig_cells_intervention$cell_type)%>%
select(patient_id, pre_post_intervention, cell_type, cell_value) %>%
group_by(cell_type) %>%
pivot_wider(names_from = pre_post_intervention,
values_from = cell_value) %>%
mutate(yellow_FC = post_Yellow/pre_Yellow,
red_FC = post_Red/pre_Red,
post_intervention_FC = post_Red/post_Yellow)
intervention_sigcells_subset %>%
summarize(mean_yellow_FC = mean(yellow_FC, na.rm = TRUE),
mean_red_FC = mean(red_FC),
mean_intervention_FC = mean(post_intervention_FC, na.rm = TRUE))
## # A tibble: 2 × 4
## cell_type mean_yellow_FC mean_red_FC mean_intervention_FC
## <chr> <dbl> <dbl> <dbl>
## 1 x38_cd16_nk_cells 1.25 1.08 0.427
## 2 x40_cd14_mdsc_mono 0.976 1.52 43.7
Very large fold change of approx 44 units higher post-Red vs. post-Yellow for MDSC cell type.
Let’s look at the correlation between significant immuno outcomes (from each comparison) and carotenoids. I also plan to correlate these with urinary soy isoflavones and their metabolites, along with other significant urinary metabolites. I will perform a log2 fold change transformation for pre-post comparison and intervention comparisons.
Pre vs post yellow
wrangling
# create df with pre-interventions data
meta_table_pre_subset <- imputed_meta_table %>%
filter(pre_post == "pre") %>%
mutate_at(11:ncol(.), log2)
# create df with post-interventions data
meta_table_post_subset <- imputed_meta_table %>%
filter(pre_post == "post") %>%
mutate_at(11:ncol(.), log2)
# subtract pre df from post df
meta_table_post_pre_differences <- meta_table_post_subset[,11:ncol(imputed_meta_table)] - meta_table_pre_subset[,11:ncol(imputed_meta_table)]
# add metadata back in and organize so that metadata is at beginning of df
meta_table_post_pre_differences <- meta_table_post_pre_differences %>%
mutate(patient_id = meta_table_post_subset$patient_id,
intervention = meta_table_post_subset$intervention,
age_at_enrollment = meta_table_post_subset$age_at_enrollment,
sex = meta_table_post_subset$sex,
bmi_at_enrollment = meta_table_post_subset$bmi_at_enrollment) %>%
relocate(patient_id, intervention, sex, age_at_enrollment, bmi_at_enrollment)
Corr table
correlation_yellow <- meta_table_post_pre_differences %>%
filter(intervention == "Yellow") %>%
filter(patient_id != 6112) %>% # remove this subj since they were an outlier in control juice for total lyc
select(total_lyc, all_of(sig_cells_yellow$cell_type), all_of(sig_cytokines_yellow$cytokine)) %>%
correlate(method = "spearman")
kable(correlation_yellow, format = "markdown", digits = 3)
| term | total_lyc | x05_cd4_cd8_cd8_t_cells | x08_cd45ro_cd45ra_em_cd8 | x09_cd45r0_cd45ra_te_cd8 | x15_cd45ro_cd45ra_te_cd4 | x26_cd27_ig_d_naive_b_cells | x37_cd56_cd161_cd123_nk_cells |
|---|---|---|---|---|---|---|---|
| total_lyc | NA | -0.627 | -0.545 | -0.573 | -0.291 | -0.473 | 0.000 |
| x05_cd4_cd8_cd8_t_cells | -0.627 | NA | 0.273 | 0.609 | 0.455 | 0.673 | -0.500 |
| x08_cd45ro_cd45ra_em_cd8 | -0.545 | 0.273 | NA | 0.518 | 0.418 | 0.100 | 0.209 |
| x09_cd45r0_cd45ra_te_cd8 | -0.573 | 0.609 | 0.518 | NA | 0.800 | 0.555 | -0.264 |
| x15_cd45ro_cd45ra_te_cd4 | -0.291 | 0.455 | 0.418 | 0.800 | NA | 0.491 | -0.009 |
| x26_cd27_ig_d_naive_b_cells | -0.473 | 0.673 | 0.100 | 0.555 | 0.491 | NA | -0.436 |
| x37_cd56_cd161_cd123_nk_cells | 0.000 | -0.500 | 0.209 | -0.264 | -0.009 | -0.436 | NA |
correlation_yellow %>%
rearrange(absolute = FALSE) %>%
shave() %>%
rplot(shape = 19,
colors = c("red", "green"),
print_cor = TRUE) +
theme(axis.text.x = element_text(angle = 90))
correlation_red <- meta_table_post_pre_differences %>%
filter(intervention == "Red") %>%
select(total_lyc, all_of(sig_cells_red$cell_type), all_of(sig_cytokines_red$cytokine)) %>%
correlate(method = "spearman")
kable(correlation_red, format = "markdown", digits = 3)
| term | total_lyc | x08_cd45ro_cd45ra_em_cd8 | x25_cd19_cd3_b_cells | x26_cd27_ig_d_naive_b_cells | gm_csf | il_12p70 | il_5 |
|---|---|---|---|---|---|---|---|
| total_lyc | NA | 0.091 | 0.140 | 0.035 | 0.007 | -0.270 | -0.189 |
| x08_cd45ro_cd45ra_em_cd8 | 0.091 | NA | 0.168 | 0.238 | 0.427 | 0.147 | 0.070 |
| x25_cd19_cd3_b_cells | 0.140 | 0.168 | NA | 0.958 | 0.042 | -0.021 | -0.161 |
| x26_cd27_ig_d_naive_b_cells | 0.035 | 0.238 | 0.958 | NA | 0.007 | -0.102 | -0.224 |
| gm_csf | 0.007 | 0.427 | 0.042 | 0.007 | NA | 0.361 | -0.028 |
| il_12p70 | -0.270 | 0.147 | -0.021 | -0.102 | 0.361 | NA | 0.666 |
| il_5 | -0.189 | 0.070 | -0.161 | -0.224 | -0.028 | 0.666 | NA |
correlation_red %>%
rearrange(absolute = FALSE) %>%
shave() %>%
rplot(shape = 19,
colors = c("red", "green"),
print_cor = TRUE) +
theme(axis.text.x = element_text(angle = 90))
wrangling
# create df with pre-interventions data
meta_table_postY_subset <- imputed_meta_table %>%
filter(pre_post_intervention == "post_Yellow") %>%
mutate_at(11:ncol(.), log2)
# create df with post-interventions data
meta_table_postR_subset <- imputed_meta_table %>%
filter(pre_post_intervention == "post_Red") %>%
mutate_at(11:ncol(.), log2)
# subtract pre df from post df
meta_table_intervention_differences <- meta_table_postR_subset[,11:ncol(imputed_meta_table)] - meta_table_postY_subset[,11:ncol(imputed_meta_table)]
# add metadata back in and organize so that metadata is at beginning of df
meta_table_intervention_differences <- meta_table_intervention_differences %>%
mutate(patient_id = meta_table_postR_subset$patient_id,
intervention = meta_table_postR_subset$intervention,
age_at_enrollment = meta_table_postR_subset$age_at_enrollment,
sex = meta_table_postR_subset$sex,
bmi_at_enrollment = meta_table_postR_subset$bmi_at_enrollment) %>%
relocate(patient_id, intervention, sex, age_at_enrollment, bmi_at_enrollment)
correlation_intervention <- meta_table_intervention_differences %>%
filter(intervention == "Red") %>%
select(total_lyc, all_of(sig_cells_intervention$cell_type), all_of(sig_cytokines_intervention$cytokine)) %>%
correlate(method = "spearman")
kable(correlation_intervention, format = "markdown", digits = 3)
| term | total_lyc | x38_cd16_nk_cells | x40_cd14_mdsc_mono |
|---|---|---|---|
| total_lyc | NA | -0.175 | -0.140 |
| x38_cd16_nk_cells | -0.175 | NA | 0.126 |
| x40_cd14_mdsc_mono | -0.140 | 0.126 | NA |
correlation_intervention %>%
rearrange(absolute = FALSE) %>%
shave() %>%
rplot(shape = 19,
colors = c("red", "green"),
print_cor = TRUE)
Let’s take all of the significant outcomes and even add all of the metadata
correlation_all_pre_post_sig <- meta_table_post_pre_differences %>%
select(2:5, total_lyc, total_carotenoids, a_carotene, b_carotene, lutein, zeaxanthin, b_cryptoxanthin, all_of(sig_cells_red$cell_type), all_of(sig_cytokines_red$cytokine), all_of(sig_cells_yellow$cell_type), all_of(sig_cytokines_yellow$cytokine), all_of(sig_cells_intervention$cell_type), all_of(sig_cytokines_intervention$cytokine)) %>%
correlate(method = "spearman")
kable(correlation_all_pre_post_sig, format = "markdown", digits = 3)
| term | age_at_enrollment | bmi_at_enrollment | total_lyc | total_carotenoids | a_carotene | b_carotene | lutein | zeaxanthin | b_cryptoxanthin | x08_cd45ro_cd45ra_em_cd8 | x25_cd19_cd3_b_cells | x26_cd27_ig_d_naive_b_cells | gm_csf | il_12p70 | il_5 | x05_cd4_cd8_cd8_t_cells | x09_cd45r0_cd45ra_te_cd8 | x15_cd45ro_cd45ra_te_cd4 | x37_cd56_cd161_cd123_nk_cells | x38_cd16_nk_cells | x40_cd14_mdsc_mono |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| age_at_enrollment | NA | 0.490 | -0.134 | -0.101 | -0.234 | -0.143 | -0.178 | 0.040 | 0.110 | 0.256 | 0.026 | 0.087 | -0.223 | -0.112 | 0.033 | 0.387 | 0.038 | 0.192 | 0.106 | 0.007 | -0.169 |
| bmi_at_enrollment | 0.490 | NA | -0.146 | -0.047 | 0.143 | -0.096 | 0.023 | 0.077 | -0.061 | 0.094 | -0.070 | 0.000 | -0.230 | 0.007 | -0.131 | 0.366 | 0.267 | 0.077 | 0.201 | -0.242 | 0.218 |
| total_lyc | -0.134 | -0.146 | NA | 0.819 | 0.120 | 0.676 | 0.161 | -0.365 | 0.082 | -0.363 | 0.023 | 0.036 | -0.157 | -0.490 | -0.424 | -0.348 | -0.262 | -0.163 | -0.267 | 0.083 | 0.029 |
| total_carotenoids | -0.101 | -0.047 | 0.819 | NA | 0.378 | 0.847 | 0.527 | -0.175 | 0.211 | -0.226 | -0.143 | -0.117 | -0.217 | -0.347 | -0.365 | -0.350 | -0.201 | -0.210 | -0.035 | 0.051 | 0.141 |
| a_carotene | -0.234 | 0.143 | 0.120 | 0.378 | NA | 0.530 | 0.437 | 0.194 | 0.117 | 0.409 | -0.195 | -0.232 | 0.218 | 0.217 | 0.220 | 0.054 | 0.354 | 0.244 | 0.290 | 0.073 | 0.369 |
| b_carotene | -0.143 | -0.096 | 0.676 | 0.847 | 0.530 | NA | 0.396 | -0.300 | -0.092 | -0.063 | -0.340 | -0.347 | -0.036 | -0.131 | -0.124 | -0.303 | -0.092 | -0.088 | 0.123 | 0.183 | 0.326 |
| lutein | -0.178 | 0.023 | 0.161 | 0.527 | 0.437 | 0.396 | NA | 0.288 | 0.150 | 0.245 | -0.092 | -0.011 | 0.063 | 0.004 | -0.205 | 0.058 | -0.042 | 0.014 | 0.089 | -0.106 | 0.004 |
| zeaxanthin | 0.040 | 0.077 | -0.365 | -0.175 | 0.194 | -0.300 | 0.288 | NA | 0.170 | 0.253 | 0.117 | 0.165 | -0.157 | 0.400 | 0.268 | 0.202 | -0.109 | 0.003 | 0.096 | -0.075 | -0.219 |
| b_cryptoxanthin | 0.110 | -0.061 | 0.082 | 0.211 | 0.117 | -0.092 | 0.150 | 0.170 | NA | 0.124 | 0.148 | 0.182 | -0.121 | -0.077 | -0.124 | -0.028 | 0.102 | -0.003 | 0.153 | 0.038 | 0.068 |
| x08_cd45ro_cd45ra_em_cd8 | 0.256 | 0.094 | -0.363 | -0.226 | 0.409 | -0.063 | 0.245 | 0.253 | 0.124 | NA | 0.101 | 0.088 | 0.312 | 0.350 | 0.411 | 0.523 | 0.441 | 0.677 | 0.378 | 0.083 | 0.167 |
| x25_cd19_cd3_b_cells | 0.026 | -0.070 | 0.023 | -0.143 | -0.195 | -0.340 | -0.092 | 0.117 | 0.148 | 0.101 | NA | 0.963 | -0.070 | -0.044 | 0.063 | 0.462 | 0.382 | 0.530 | -0.323 | -0.183 | -0.213 |
| x26_cd27_ig_d_naive_b_cells | 0.087 | 0.000 | 0.036 | -0.117 | -0.232 | -0.347 | -0.011 | 0.165 | 0.182 | 0.088 | 0.963 | NA | -0.047 | -0.073 | -0.073 | 0.581 | 0.308 | 0.512 | -0.406 | -0.288 | -0.324 |
| gm_csf | -0.223 | -0.230 | -0.157 | -0.217 | 0.218 | -0.036 | 0.063 | -0.157 | -0.121 | 0.312 | -0.070 | -0.047 | NA | 0.267 | 0.136 | 0.215 | -0.002 | 0.244 | 0.023 | -0.046 | -0.098 |
| il_12p70 | -0.112 | 0.007 | -0.490 | -0.347 | 0.217 | -0.131 | 0.004 | 0.400 | -0.077 | 0.350 | -0.044 | -0.073 | 0.267 | NA | 0.556 | 0.114 | 0.204 | 0.125 | 0.584 | 0.341 | 0.142 |
| il_5 | 0.033 | -0.131 | -0.424 | -0.365 | 0.220 | -0.124 | -0.205 | 0.268 | -0.124 | 0.411 | 0.063 | -0.073 | 0.136 | 0.556 | NA | 0.028 | 0.320 | 0.273 | 0.402 | 0.568 | 0.186 |
| x05_cd4_cd8_cd8_t_cells | 0.387 | 0.366 | -0.348 | -0.350 | 0.054 | -0.303 | 0.058 | 0.202 | -0.028 | 0.523 | 0.462 | 0.581 | 0.215 | 0.114 | 0.028 | NA | 0.461 | 0.678 | -0.158 | -0.428 | -0.154 |
| x09_cd45r0_cd45ra_te_cd8 | 0.038 | 0.267 | -0.262 | -0.201 | 0.354 | -0.092 | -0.042 | -0.109 | 0.102 | 0.441 | 0.382 | 0.308 | -0.002 | 0.204 | 0.320 | 0.461 | NA | 0.566 | 0.234 | 0.015 | 0.258 |
| x15_cd45ro_cd45ra_te_cd4 | 0.192 | 0.077 | -0.163 | -0.210 | 0.244 | -0.088 | 0.014 | 0.003 | -0.003 | 0.677 | 0.530 | 0.512 | 0.244 | 0.125 | 0.273 | 0.678 | 0.566 | NA | 0.110 | 0.041 | -0.012 |
| x37_cd56_cd161_cd123_nk_cells | 0.106 | 0.201 | -0.267 | -0.035 | 0.290 | 0.123 | 0.089 | 0.096 | 0.153 | 0.378 | -0.323 | -0.406 | 0.023 | 0.584 | 0.402 | -0.158 | 0.234 | 0.110 | NA | 0.589 | 0.464 |
| x38_cd16_nk_cells | 0.007 | -0.242 | 0.083 | 0.051 | 0.073 | 0.183 | -0.106 | -0.075 | 0.038 | 0.083 | -0.183 | -0.288 | -0.046 | 0.341 | 0.568 | -0.428 | 0.015 | 0.041 | 0.589 | NA | 0.103 |
| x40_cd14_mdsc_mono | -0.169 | 0.218 | 0.029 | 0.141 | 0.369 | 0.326 | 0.004 | -0.219 | 0.068 | 0.167 | -0.213 | -0.324 | -0.098 | 0.142 | 0.186 | -0.154 | 0.258 | -0.012 | 0.464 | 0.103 | NA |
correlation_all_pre_post_sig %>%
rearrange(absolute = FALSE) %>%
shave() %>%
rplot(shape = 19,
print_cor = TRUE,
colors = c("red", "green")) +
theme(axis.text.x = element_text(angle = 90))
correlation_all_intervention_sig <- meta_table_intervention_differences %>%
select(2:5, total_lyc, total_carotenoids, a_carotene, b_carotene, lutein, zeaxanthin, b_cryptoxanthin, all_of(sig_cells_red$cell_type), all_of(sig_cytokines_red$cytokine), all_of(sig_cells_yellow$cell_type), all_of(sig_cytokines_yellow$cytokine), all_of(sig_cells_intervention$cell_type), all_of(sig_cytokines_intervention$cytokine)) %>%
correlate(method = "spearman")
kable(correlation_all_intervention_sig, format = "markdown", digits = 3)
| term | age_at_enrollment | bmi_at_enrollment | total_lyc | total_carotenoids | a_carotene | b_carotene | lutein | zeaxanthin | b_cryptoxanthin | x08_cd45ro_cd45ra_em_cd8 | x25_cd19_cd3_b_cells | x26_cd27_ig_d_naive_b_cells | gm_csf | il_12p70 | il_5 | x05_cd4_cd8_cd8_t_cells | x09_cd45r0_cd45ra_te_cd8 | x15_cd45ro_cd45ra_te_cd4 | x37_cd56_cd161_cd123_nk_cells | x38_cd16_nk_cells | x40_cd14_mdsc_mono |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| age_at_enrollment | NA | 0.490 | 0.056 | -0.007 | -0.238 | 0.343 | 0.014 | 0.455 | -0.203 | 0.105 | 0.538 | 0.510 | 0.322 | 0.385 | 0.224 | 0.035 | 0.266 | 0.294 | -0.147 | -0.105 | -0.147 |
| bmi_at_enrollment | 0.490 | NA | -0.126 | -0.098 | -0.049 | 0.189 | -0.231 | 0.119 | -0.119 | -0.406 | 0.063 | 0.021 | 0.287 | 0.252 | 0.028 | -0.490 | 0.161 | -0.357 | 0.615 | 0.147 | 0.350 |
| total_lyc | 0.056 | -0.126 | NA | 0.797 | 0.545 | 0.434 | 0.741 | 0.643 | 0.587 | -0.294 | -0.280 | -0.196 | 0.301 | 0.091 | 0.329 | -0.294 | -0.273 | 0.280 | -0.119 | -0.175 | -0.140 |
| total_carotenoids | -0.007 | -0.098 | 0.797 | NA | 0.601 | 0.671 | 0.944 | 0.587 | 0.867 | -0.385 | -0.161 | -0.203 | 0.182 | 0.070 | 0.413 | -0.510 | -0.343 | -0.028 | 0.056 | -0.294 | -0.063 |
| a_carotene | -0.238 | -0.049 | 0.545 | 0.601 | NA | 0.490 | 0.434 | 0.098 | 0.727 | -0.441 | -0.427 | -0.399 | -0.161 | -0.175 | 0.133 | -0.294 | -0.888 | -0.294 | 0.182 | 0.049 | -0.098 |
| b_carotene | 0.343 | 0.189 | 0.434 | 0.671 | 0.490 | NA | 0.608 | 0.245 | 0.497 | -0.056 | -0.105 | -0.196 | -0.126 | -0.056 | 0.552 | -0.322 | -0.350 | -0.021 | -0.028 | -0.469 | -0.287 |
| lutein | 0.014 | -0.231 | 0.741 | 0.944 | 0.434 | 0.608 | NA | 0.650 | 0.804 | -0.224 | -0.049 | -0.042 | 0.287 | 0.126 | 0.476 | -0.413 | -0.196 | 0.168 | -0.049 | -0.385 | -0.273 |
| zeaxanthin | 0.455 | 0.119 | 0.643 | 0.587 | 0.098 | 0.245 | 0.650 | NA | 0.399 | -0.231 | 0.231 | 0.266 | 0.629 | 0.538 | 0.329 | -0.224 | 0.140 | 0.371 | -0.028 | -0.210 | -0.035 |
| b_cryptoxanthin | -0.203 | -0.119 | 0.587 | 0.867 | 0.727 | 0.497 | 0.804 | 0.399 | NA | -0.587 | -0.175 | -0.168 | 0.070 | -0.028 | 0.091 | -0.559 | -0.615 | -0.343 | 0.182 | -0.091 | -0.035 |
| x08_cd45ro_cd45ra_em_cd8 | 0.105 | -0.406 | -0.294 | -0.385 | -0.441 | -0.056 | -0.224 | -0.231 | -0.587 | NA | -0.091 | -0.105 | -0.392 | -0.385 | 0.119 | 0.839 | 0.413 | 0.685 | -0.636 | -0.573 | -0.301 |
| x25_cd19_cd3_b_cells | 0.538 | 0.063 | -0.280 | -0.161 | -0.427 | -0.105 | -0.049 | 0.231 | -0.175 | -0.091 | NA | 0.944 | 0.503 | 0.748 | 0.175 | -0.091 | 0.273 | 0.028 | 0.049 | 0.392 | -0.252 |
| x26_cd27_ig_d_naive_b_cells | 0.510 | 0.021 | -0.196 | -0.203 | -0.399 | -0.196 | -0.042 | 0.266 | -0.168 | -0.105 | 0.944 | NA | 0.615 | 0.727 | 0.126 | -0.077 | 0.238 | 0.140 | 0.014 | 0.441 | -0.420 |
| gm_csf | 0.322 | 0.287 | 0.301 | 0.182 | -0.161 | -0.126 | 0.287 | 0.629 | 0.070 | -0.392 | 0.503 | 0.615 | NA | 0.797 | 0.420 | -0.490 | 0.357 | 0.245 | 0.350 | 0.392 | -0.189 |
| il_12p70 | 0.385 | 0.252 | 0.091 | 0.070 | -0.175 | -0.056 | 0.126 | 0.538 | -0.028 | -0.385 | 0.748 | 0.727 | 0.797 | NA | 0.420 | -0.371 | 0.217 | 0.049 | 0.343 | 0.448 | -0.133 |
| il_5 | 0.224 | 0.028 | 0.329 | 0.413 | 0.133 | 0.552 | 0.476 | 0.329 | 0.091 | 0.119 | 0.175 | 0.126 | 0.420 | 0.420 | NA | -0.259 | 0.140 | 0.434 | 0.000 | -0.112 | -0.413 |
| x05_cd4_cd8_cd8_t_cells | 0.035 | -0.490 | -0.294 | -0.510 | -0.294 | -0.322 | -0.413 | -0.224 | -0.559 | 0.839 | -0.091 | -0.077 | -0.490 | -0.371 | -0.259 | NA | 0.147 | 0.531 | -0.657 | -0.357 | -0.147 |
| x09_cd45r0_cd45ra_te_cd8 | 0.266 | 0.161 | -0.273 | -0.343 | -0.888 | -0.350 | -0.196 | 0.140 | -0.615 | 0.413 | 0.273 | 0.238 | 0.357 | 0.217 | 0.140 | 0.147 | NA | 0.455 | -0.084 | -0.161 | 0.168 |
| x15_cd45ro_cd45ra_te_cd4 | 0.294 | -0.357 | 0.280 | -0.028 | -0.294 | -0.021 | 0.168 | 0.371 | -0.343 | 0.685 | 0.028 | 0.140 | 0.245 | 0.049 | 0.434 | 0.531 | 0.455 | NA | -0.643 | -0.371 | -0.434 |
| x37_cd56_cd161_cd123_nk_cells | -0.147 | 0.615 | -0.119 | 0.056 | 0.182 | -0.028 | -0.049 | -0.028 | 0.182 | -0.636 | 0.049 | 0.014 | 0.350 | 0.343 | 0.000 | -0.657 | -0.084 | -0.643 | NA | 0.350 | 0.217 |
| x38_cd16_nk_cells | -0.105 | 0.147 | -0.175 | -0.294 | 0.049 | -0.469 | -0.385 | -0.210 | -0.091 | -0.573 | 0.392 | 0.441 | 0.392 | 0.448 | -0.112 | -0.357 | -0.161 | -0.371 | 0.350 | NA | 0.126 |
| x40_cd14_mdsc_mono | -0.147 | 0.350 | -0.140 | -0.063 | -0.098 | -0.287 | -0.273 | -0.035 | -0.035 | -0.301 | -0.252 | -0.420 | -0.189 | -0.133 | -0.413 | -0.147 | 0.168 | -0.434 | 0.217 | 0.126 | NA |
correlation_all_intervention_sig %>%
rearrange(absolute = FALSE) %>%
shave() %>%
rplot(shape = 19,
print_cor = TRUE,
colors = c("red", "green")) +
theme(axis.text.x = element_text(angle = 90))
Here, I am using my raw data and allowing the pheatmap package to perform z-scaling.
Z-scoring standardizes the data by this calculation: (individual value within outcome - mean of outcome) / (std dev). The pheatmap package does this automatically with the call scale = “row” or “column”
cytokine_heatmap_data <- imputed_meta_table %>%
unite("subject_pre_post_intervention", patient_id, pre_post_intervention, sep = "_", remove = FALSE) %>%
dplyr::select(patient_id, subject_pre_post_intervention, pre_post, all_of(sig_cytokines_red$cytokine), all_of(sig_cytokines_yellow$cytokine), all_of(sig_cytokines_intervention$cytokine))
cytokine_heatmap <-
pheatmap(cytokine_heatmap_data[,-c(1:3)],
scale = "column", # z-scaling
cluster_rows = TRUE,
cutree_rows = 6,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D2",
labels_row = cytokine_heatmap_data$subject_pre_post_intervention,
color = colorRampPalette(c("#67a9cf", "#f7f7f7", "#ef8a62"))(16),
main = "Heatmap of significant immuno outcomes across timepoints \nby paired T-tests \nBenjamoni-Hochberg corrected p-values > 0.05 \nC18 (-)")
cytokines_all_heatmap_data <- imputed_meta_table %>%
unite("subject_pre_post_intervention", patient_id, pre_post_intervention, sep = "_", remove = FALSE) %>%
dplyr::select(patient_id, subject_pre_post_intervention, pre_post, all_of(cytokines_long$cytokine))
cytokines_all_heatmap <-
pheatmap(cytokines_all_heatmap_data[,-c(1:3)],
scale = "column", # z-scaling
cluster_rows = TRUE,
cutree_rows = 4,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D",
labels_row = cytokines_all_heatmap_data$subject_pre_post_intervention,
color = colorRampPalette(c("#67a9cf", "#f7f7f7", "#ef8a62"))(16),
main = "Heatmap of all cytokines from multiplex assay")
I want to try a different technique here to see if we can capture variation by bringing the values to a visually meaningful range by normalizing values against “controls”. Each subject serves as their own control pre-intervention (and also when they are on the low carotenoid/soy intervention (Yellow)). These dfs will be different from the log2 fold change dfs calculated for correlation plots. In this case, we will take every treatment and divide it by its control (so control should be end up being 0)
# make a string of metadata and of immuno outcomes for ease
str_meta <- colnames(imputed_meta_table[,1:10])
str_immuno <- colnames(imputed_meta_table[,11:82])
# Extract pre-intervention data and set immuno outcome columns to 0
pre_data_zeros <- imputed_meta_table %>%
filter(pre_post == "pre") %>%
mutate(across(all_of(str_immuno), ~ 0))
# Combine pre-intervention zeros with log fold change table
normalized_to_pre <- meta_table_post_pre_differences %>%
mutate(pre_post = "normalized_POST") %>%
bind_rows(pre_data_zeros %>%
select(str_meta[c(1:4, 6)], all_of(str_immuno)) %>%
mutate(pre_post = "normalized_PRE")) %>%
unite("pre_post_intervention", pre_post, intervention, sep = "_", remove = FALSE)
cytokines_normlog2FC_heatmap_data <- normalized_to_pre %>%
unite("subj_pre_post_intervention", patient_id, pre_post_intervention, sep = "_", remove = FALSE) %>%
filter(pre_post == "normalized_POST") %>%
dplyr::select(subj_pre_post_intervention, all_of(cytokines_long$cytokine))
cytokines_normtopre_heatmap <-
pheatmap(cytokines_normlog2FC_heatmap_data[,-1],
cluster_rows = TRUE,
cutree_rows = 2,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D",
labels_row = cytokines_normlog2FC_heatmap_data$subj_pre_post_intervention,
color = colorRampPalette(c("#67a9cf", "#f7f7f7", "#ef8a62"))(16),
main = "Heatmap of all cytokines from multiplex assay")
cytokines2_normlog2FC_heatmap_data <- normalized_to_pre %>%
unite("subj_pre_post_intervention", patient_id, pre_post_intervention, sep = "_", remove = FALSE) %>%
dplyr::select(subj_pre_post_intervention, all_of(cytokines_long$cytokine))
cytokines2_normtopre_heatmap <-
pheatmap(cytokines2_normlog2FC_heatmap_data[,-1],
cluster_rows = TRUE,
cutree_rows = 2,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D",
labels_row = cytokines2_normlog2FC_heatmap_data$subj_pre_post_intervention,
color = colorRampPalette(c("#67a9cf", "#f7f7f7", "#ef8a62"))(16),
main = "Heatmap of all cytokines from multiplex assay")
yellow_norm_to_pre <- normalized_to_pre %>%
filter(intervention == "Yellow")
cytokines_normlog2FC_yellow_heatmap_data <- yellow_norm_to_pre %>%
unite("subject_pre_post", patient_id, pre_post, sep = "_", remove = FALSE) %>%
dplyr::select(subject_pre_post, all_of(cytokines_long$cytokine))
cytokines_yellow_heatmap <-
pheatmap(cytokines_normlog2FC_yellow_heatmap_data[,-1],
cluster_rows = TRUE,
cutree_rows = 2,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D",
labels_row = cytokines_normlog2FC_yellow_heatmap_data$subject_pre_post,
color = colorRampPalette(c("#67a9cf", "#f7f7f7", "#ef8a62"))(16),
main = "Heatmap of all cytokines from multiplex assay")
red_norm_to_pre <- normalized_to_pre %>%
filter(intervention == "Red")
cytokines_normlog2FC_red_heatmap_data <- red_norm_to_pre %>%
unite("subject_pre_post", patient_id, pre_post, sep = "_", remove = FALSE) %>%
dplyr::select(subject_pre_post, all_of(cytokines_long$cytokine))
cytokines_red_heatmap <-
pheatmap(cytokines_normlog2FC_red_heatmap_data[,-1],
cluster_rows = TRUE,
cutree_rows = 3,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D",
labels_row = cytokines_normlog2FC_red_heatmap_data$subject_pre_post,
color = colorRampPalette(c("#67a9cf", "#f7f7f7", "#ef8a62"))(16),
main = "Heatmap of all cytokines from multiplex assay")
cell_heatmap_data <- imputed_meta_table %>%
unite("subject_pre_post_intervention", patient_id, pre_post_intervention, sep = "_", remove = FALSE) %>%
dplyr::select(patient_id, subject_pre_post_intervention, pre_post, all_of(sig_cells_red$cell_type), all_of(sig_cells_yellow$cell_type), all_of(sig_cells_intervention$cell_type))
cell_heatmap <-
pheatmap(cell_heatmap_data[,-c(1:3)],
scale = "column",
cluster_rows = TRUE,
cutree_rows = 6,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "ward.D2",
labels_row = cell_heatmap_data$subject_pre_post_intervention,
color = colorRampPalette(c("#67a9cf", "#f7f7f7", "#ef8a62"))(16),
main = "Heatmap of significant immuno outcomes across timepoints \nby paired T-tests \nBenjamoni-Hochberg corrected p-values > 0.05 \nC18 (-)")
There are no sequence effects for any of the outcomes
Carotenoids
Cytokines
Immune cell types